Skip to content

Commit 8a81ecb

Browse files
slipchenkoslipchenko
committed
working and clean version before checking with qchem
1 parent 289ff8f commit 8a81ecb

7 files changed

Lines changed: 66104 additions & 399 deletions

File tree

‎fraglib/h2o_reduced.efp‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
RUNTYP=MAKEFP EFFECTIVE FRAGMENT POTENTIAL DATA FOLLOWS...
2+
FRAGNAMEEFP GENERATED AT Tue Mar 27 16:14:27 2012
3+
$H2O_REDUCED_L
4+
EFP DATA FOR FRAGNAME SCFTYP=RHF ... GENERATED WITH BASIS SET=6-311++G(3df,2p)
5+
COORDINATES (BOHR)
6+
A01O1 0.0000000000 0.1191094785 0.0000000000 15.9949100 8.0
7+
A02H2 -1.4223059670 -0.9451766865 0.0000000000 1.0078250 1.0
8+
A03H3 1.4223059670 -0.9451766865 0.0000000000 1.0078250 1.0
9+
BO21 -0.7111529835 -0.4130336040 0.0000000000 0.0000000 0.0
10+
BO31 0.7111529835 -0.4130336040 0.0000000000 0.0000000 0.0
11+
STOP
12+
MONOPOLES
13+
A01O1 -8.4173680991 8.00000
14+
A02H2 -0.3512187169 1.00000
15+
A03H3 -0.3512187169 1.00000
16+
BO21 -0.4400972335 0.00000
17+
BO31 -0.4400972335 0.00000
18+
STOP
19+
DIPOLES
20+
A01O1 0.0000000000 -0.4127062480 0.0000000000
21+
A02H2 -0.0027268078 -0.0014022974 0.0000000000
22+
A03H3 0.0027268078 -0.0014022974 0.0000000000
23+
BO21 0.2657818726 0.2358863756 0.0000000000
24+
BO31 -0.2657818726 0.2358863756 0.0000000000
25+
STOP
26+
POLARIZABLE POINTS
27+
CT1 -0.7491370634 -0.4890369778 0.0000000020
28+
2.6629294234 1.7670590537 0.7812922107 1.2832534969 >
29+
-0.0000000037 -0.0000000033 1.1482095947 -0.0000000144 >
30+
-0.0000000132
31+
CT2 0.7491370628 -0.4890369781 0.0000000040
32+
2.6629294254 1.7670590567 0.7812922113 -1.2832535009 >
33+
0.0000000076 -0.0000000067 -1.1482095971 0.0000000195 >
34+
-0.0000000177
35+
CT3 0.0000000011 0.3892153728 -0.4957986131
36+
1.4502428775 1.8574126776 2.7113709049 0.0000000046 >
37+
-0.0000000029 -0.9972382204 0.0000000048 -0.0000000034 >
38+
-1.0913249895
39+
CT4 -0.0000000005 0.3892153799 0.4957986071
40+
1.4502428978 1.8574127185 2.7113708742 -0.0000000023 >
41+
-0.0000000009 0.9972382303 -0.0000000024 -0.0000000017 >
42+
1.0913250020
43+
STOP
44+
$END

‎src/efp.c‎

Lines changed: 97 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -409,25 +409,35 @@ check_frag_params(const struct efp_opts *opts, struct frag *frag)
409409
{
410410
if ((opts->terms & EFP_TERM_ELEC) || (opts->terms & EFP_TERM_AI_ELEC)) {
411411
if (!frag->multipole_pts) {
412-
efp_log("electrostatic parameters are missing");
413-
return EFP_RESULT_FATAL;
412+
printf("WARNING! Multipole parameters are missing on fragment %s\n", frag->name);
413+
// efp_log("electrostatic parameters are missing");
414+
// return EFP_RESULT_FATAL;
414415
}
415416
}
416417
if ((opts->terms & EFP_TERM_POL) || (opts->terms & EFP_TERM_AI_POL)) {
417-
if (!frag->polarizable_pts || !frag->multipole_pts) {
418-
efp_log("polarization parameters are missing");
419-
return EFP_RESULT_FATAL;
418+
if (!frag->multipole_pts)
419+
printf("WARNING! Multipole parameters are missing on fragment %s\n", frag->name);
420+
if (!frag->polarizable_pts) {
421+
printf("WARNING! Polarizability parameters are missing on fragment %s\n", frag->name);
422+
423+
// efp_log("polarization parameters are missing");
424+
// return EFP_RESULT_FATAL;
420425
}
421426
}
422427
if ((opts->terms & EFP_TERM_DISP) || (opts->terms & EFP_TERM_AI_DISP)) {
423428
if (frag->dynamic_polarizable_pts == NULL) {
424-
efp_log("dispersion parameters are missing");
425-
return EFP_RESULT_FATAL;
429+
printf("WARNING! Dynamic polarizability parameters are "
430+
"missing on fragment %s\n", frag->name);
431+
432+
// efp_log("dispersion parameters are missing");
433+
// return EFP_RESULT_FATAL;
426434
}
427435
if (opts->disp_damp == EFP_DISP_DAMP_OVERLAP &&
428436
frag->n_lmo != frag->n_dynamic_polarizable_pts) {
429-
efp_log("number of polarization points does not "
430-
"match number of LMOs");
437+
//efp_log("number of polarization points does not "
438+
// "match number of LMOs");
439+
printf("FATAL ERROR! The number of polarization points does not "
440+
"match the number of LMOs of fragment %s\n", frag->name);
431441
return EFP_RESULT_FATAL;
432442
}
433443
}
@@ -436,8 +446,10 @@ check_frag_params(const struct efp_opts *opts, struct frag *frag)
436446
!frag->xr_fock_mat ||
437447
!frag->xr_wf ||
438448
!frag->lmo_centroids) {
439-
efp_log("exchange repulsion parameters are missing");
440-
return EFP_RESULT_FATAL;
449+
printf("WARNING! Exchange-repulsion parameters are "
450+
"missing on fragment %s\n", frag->name);
451+
//efp_log("exchange repulsion parameters are missing");
452+
//return EFP_RESULT_FATAL;
441453
}
442454
}
443455
return EFP_RESULT_SUCCESS;
@@ -507,30 +519,34 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to,
507519
size_t n_lmo_ij = efp->frags[i].n_lmo *
508520
efp->frags[fr_j].n_lmo;
509521

510-
s = (double *)calloc(n_lmo_ij, sizeof(double));
511-
ds = (six_t *)calloc(n_lmo_ij, sizeof(six_t));
512-
513-
if (do_xr(&efp->opts)) {
514-
double exr, ecp;
515-
516-
efp_frag_frag_xr(efp, i, fr_j,
517-
s, ds, &exr, &ecp);
518-
e_xr += exr;
519-
e_cp += ecp;
520-
521-
/* */
522-
if (efp->opts.enable_pairwise) {
523-
if (i == efp->ligand_index) {
524-
efp->pair_energies[fr_j].exchange_repulsion = exr;
525-
efp->pair_energies[fr_j].charge_penetration = ecp;
526-
}
527-
if (fr_j == efp->ligand_index) {
528-
efp->pair_energies[i].exchange_repulsion = exr;
529-
efp->pair_energies[i].charge_penetration = ecp;
522+
// ugly but simple check when not to compute exchange-repulsion and overlap screening
523+
if (n_lmo_ij > 0) {
524+
s = (double *) calloc(n_lmo_ij, sizeof(double));
525+
ds = (six_t *) calloc(n_lmo_ij, sizeof(six_t));
526+
527+
if (do_xr(&efp->opts)) {
528+
double exr, ecp;
529+
530+
efp_frag_frag_xr(efp, i, fr_j,
531+
s, ds, &exr, &ecp);
532+
e_xr += exr;
533+
e_cp += ecp;
534+
535+
/* */
536+
if (efp->opts.enable_pairwise) {
537+
if (i == efp->ligand_index) {
538+
efp->pair_energies[fr_j].exchange_repulsion = exr;
539+
efp->pair_energies[fr_j].charge_penetration = ecp;
540+
}
541+
if (fr_j == efp->ligand_index) {
542+
efp->pair_energies[i].exchange_repulsion = exr;
543+
efp->pair_energies[i].charge_penetration = ecp;
544+
}
530545
}
531546
}
532-
}
533-
if (do_elec(&efp->opts)) {
547+
}
548+
if (do_elec(&efp->opts) && efp->frags[i].n_multipole_pts > 0 &&
549+
efp->frags[fr_j].n_multipole_pts > 0) {
534550
e_elec_tmp = efp_frag_frag_elec(efp,
535551
i, fr_j);
536552
e_elec += e_elec_tmp;
@@ -542,7 +558,8 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to,
542558
efp->pair_energies[i].electrostatic = e_elec_tmp;
543559
}
544560
}
545-
if (do_disp(&efp->opts)) {
561+
if (do_disp(&efp->opts) && efp->frags[i].n_dynamic_polarizable_pts > 0 &&
562+
efp->frags[fr_j].n_dynamic_polarizable_pts > 0) {
546563
e_disp_tmp = efp_frag_frag_disp(efp,
547564
i, fr_j, s, ds);
548565
e_disp += e_disp_tmp;
@@ -554,8 +571,10 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to,
554571
efp->pair_energies[i].dispersion = e_disp_tmp;
555572
}
556573
}
557-
free(s);
558-
free(ds);
574+
if (n_lmo_ij > 0) {
575+
free(s);
576+
free(ds);
577+
}
559578
}
560579
}
561580
}
@@ -593,33 +612,35 @@ compute_two_body_crystal(struct efp *efp)
593612
six_t *ds;
594613
size_t n_lmo_ij = efp->frags[i].n_lmo *
595614
efp->frags[fr_j].n_lmo;
596-
597-
s = (double *)calloc(n_lmo_ij, sizeof(double));
598-
ds = (six_t *)calloc(n_lmo_ij, sizeof(six_t));
599-
600-
if (do_xr(&efp->opts)) {
601-
double exr, ecp;
602-
603-
efp_frag_frag_xr(efp, i, fr_j,
604-
s, ds, &exr, &ecp);
605-
e_xr += exr * factor;
606-
e_cp += ecp * factor;
607-
608-
/* */
609-
if (efp->opts.enable_pairwise) {
610-
if (i == efp->opts.ligand) {
611-
efp->pair_energies[fr_j].exchange_repulsion = exr;
612-
efp->pair_energies[fr_j].charge_penetration = ecp;
613-
}
614-
if (fr_j == efp->opts.ligand) {
615-
efp->pair_energies[i].exchange_repulsion = exr;
616-
efp->pair_energies[i].charge_penetration = ecp;
615+
// if need to compute exrep and overlap screening for this pair
616+
if (n_lmo_ij > 0) {
617+
s = (double *) calloc(n_lmo_ij, sizeof(double));
618+
ds = (six_t *) calloc(n_lmo_ij, sizeof(six_t));
619+
620+
if (do_xr(&efp->opts)) {
621+
double exr, ecp;
622+
623+
efp_frag_frag_xr(efp, i, fr_j,
624+
s, ds, &exr, &ecp);
625+
e_xr += exr * factor;
626+
e_cp += ecp * factor;
627+
628+
/* */
629+
if (efp->opts.enable_pairwise) {
630+
if (i == efp->opts.ligand) {
631+
efp->pair_energies[fr_j].exchange_repulsion = exr;
632+
efp->pair_energies[fr_j].charge_penetration = ecp;
633+
}
634+
if (fr_j == efp->opts.ligand) {
635+
efp->pair_energies[i].exchange_repulsion = exr;
636+
efp->pair_energies[i].charge_penetration = ecp;
637+
}
617638
}
618639
}
619640
}
620-
if (do_elec(&efp->opts)) {
621-
e_elec_tmp = efp_frag_frag_elec(efp,
622-
i, fr_j);
641+
if (do_elec(&efp->opts) && efp->frags[i].n_multipole_pts > 0 &&
642+
efp->frags[fr_j].n_multipole_pts > 0) {
643+
e_elec_tmp = efp_frag_frag_elec(efp, i, fr_j);
623644
e_elec += e_elec_tmp * factor;
624645

625646
/* */
@@ -630,9 +651,9 @@ compute_two_body_crystal(struct efp *efp)
630651
efp->pair_energies[i].electrostatic = e_elec_tmp;
631652
}
632653
}
633-
if (do_disp(&efp->opts)) {
634-
e_disp_tmp = efp_frag_frag_disp(efp,
635-
i, fr_j, s, ds);
654+
if (do_disp(&efp->opts) && efp->frags[i].n_dynamic_polarizable_pts > 0 &&
655+
efp->frags[fr_j].n_dynamic_polarizable_pts > 0) {
656+
e_disp_tmp = efp_frag_frag_disp(efp, i, fr_j, s, ds);
636657
e_disp += e_disp_tmp * factor;
637658
/* */
638659
if (efp->opts.enable_pairwise) {
@@ -642,9 +663,10 @@ compute_two_body_crystal(struct efp *efp)
642663
efp->pair_energies[i].dispersion = e_disp_tmp;
643664
}
644665
}
645-
free(s);
646-
free(ds);
647-
666+
if (n_lmo_ij > 0) {
667+
free(s);
668+
free(ds);
669+
}
648670
}
649671
}
650672
}
@@ -1284,17 +1306,20 @@ efp_compute(struct efp *efp, int do_gradient)
12841306

12851307
assert(efp);
12861308

1309+
static int efp_counter = 0;
1310+
12871311
if (efp->grad == NULL) {
12881312
efp_log("call efp_prepare after all fragments are added");
12891313
return EFP_RESULT_FATAL;
12901314
}
12911315

12921316
efp->do_gradient = do_gradient;
12931317

1294-
if ((res = check_params(efp))) {
1295-
efp_log("check_params() failure");
1296-
return res;
1297-
}
1318+
if (efp_counter == 0)
1319+
if ((res = check_params(efp))) {
1320+
efp_log("check_params() failure");
1321+
return res;
1322+
}
12981323

12991324
memset(&efp->energy, 0, sizeof(efp->energy));
13001325
memset(&efp->stress, 0, sizeof(efp->stress));
@@ -1345,6 +1370,8 @@ efp_compute(struct efp *efp, int do_gradient)
13451370
efp->energy.ai_dispersion +
13461371
efp->energy.exchange_repulsion;
13471372

1373+
efp_counter++;
1374+
13481375
return EFP_RESULT_SUCCESS;
13491376
}
13501377

0 commit comments

Comments
 (0)