1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
50 #include "nonbonded.h"
65 #include "mpelogging.h"
80 gmx_grppairener_t
*grppener
,
82 gmx_bool bDoLongRange
,
89 GMX_MPE_LOG(ev_ns_start
);
90 if (!fr
->ns
.nblist_initialized
)
92 init_neighbor_list(fp
, fr
, md
->homenr
);
98 nsearch
= search_neighbours(fp
,fr
,x
,box
,top
,groups
,cr
,nrnb
,md
,
99 lambda
,dvdlambda
,grppener
,
100 bFillGrid
,bDoLongRange
,
103 fprintf(debug
,"nsearch = %d\n",nsearch
);
105 /* Check whether we have to do dynamic load balancing */
106 /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
107 count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
108 &(top->idef),opts->ngener);
110 if (fr
->ns
.dump_nl
> 0)
111 dump_nblist(fp
,cr
,fr
,fr
->ns
.dump_nl
);
113 GMX_MPE_LOG(ev_ns_finish
);
116 void do_force_lowlevel(FILE *fplog
, gmx_large_int_t step
,
117 t_forcerec
*fr
, t_inputrec
*ir
,
118 t_idef
*idef
, t_commrec
*cr
,
119 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
122 rvec x
[], history_t
*hist
,
124 gmx_enerdata_t
*enerd
,
141 gmx_bool bDoEpot
,bSepDVDL
,bSB
;
145 real dvdlambda
,Vsr
,Vlr
,Vcorr
=0,vdip
,vcharge
;
149 gmx_enerdata_t ed_lam
;
154 double t0
=0.0,t1
,t2
,t3
; /* time measurement for coarse load balancing */
157 #define PRINT_SEPDVDL(s,v,dvdl) if (bSepDVDL) fprintf(fplog,sepdvdlformat,s,v,dvdl);
159 GMX_MPE_LOG(ev_force_start
);
160 set_pbc(&pbc
,fr
->ePBC
,box
);
163 for(i
=0; (i
<DIM
); i
++)
165 box_size
[i
]=box
[i
][i
];
168 bSepDVDL
=(fr
->bSepDVDL
&& do_per_step(step
,ir
->nstlog
));
171 /* do QMMM first if requested */
174 enerd
->term
[F_EQM
] = calculate_QMMM(cr
,x
,f
,fr
,md
);
179 fprintf(fplog
,"Step %s: non-bonded V and dVdl for node %d:\n",
180 gmx_step_str(step
,buf
),cr
->nodeid
);
183 /* Call the short range functions all in one go. */
184 GMX_MPE_LOG(ev_do_fnbf_start
);
189 /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
190 #define TAKETIME FALSE
193 MPI_Barrier(cr
->mpi_comm_mygroup
);
200 dvdlambda
= do_walls(ir
,fr
,box
,md
,x
,f
,lambda
,
201 enerd
->grpp
.ener
[egLJSR
],nrnb
);
202 PRINT_SEPDVDL("Walls",0.0,dvdlambda
);
203 enerd
->dvdl_lin
+= dvdlambda
;
206 /* If doing GB, reset dvda and calculate the Born radii */
207 if (ir
->implicit_solvent
)
209 /* wallcycle_start(wcycle,ewcGB); */
211 for(i
=0;i
<born
->nr
;i
++)
218 calc_gb_rad(cr
,fr
,ir
,top
,atype
,x
,&(fr
->gblist
),born
,md
,nrnb
);
221 /* wallcycle_stop(wcycle, ewcGB); */
226 if (flags
& GMX_FORCE_FORCES
)
228 donb_flags
|= GMX_DONB_FORCES
;
230 do_nonbonded(cr
,fr
,x
,f
,md
,excl
,
232 enerd
->grpp
.ener
[egBHAMSR
] :
233 enerd
->grpp
.ener
[egLJSR
],
234 enerd
->grpp
.ener
[egCOULSR
],
235 enerd
->grpp
.ener
[egGB
],box_size
,nrnb
,
236 lambda
,&dvdlambda
,-1,-1,donb_flags
);
237 /* If we do foreign lambda and we have soft-core interactions
238 * we have to recalculate the (non-linear) energies contributions.
240 if (ir
->n_flambda
> 0 && (flags
& GMX_FORCE_DHDL
) && ir
->sc_alpha
!= 0)
242 init_enerdata(mtop
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,&ed_lam
);
244 for(i
=0; i
<enerd
->n_lambda
; i
++)
246 lam_i
= (i
==0 ? lambda
: ir
->flambda
[i
-1]);
248 reset_enerdata(&ir
->opts
,fr
,TRUE
,&ed_lam
,FALSE
);
249 do_nonbonded(cr
,fr
,x
,f
,md
,excl
,
251 ed_lam
.grpp
.ener
[egBHAMSR
] :
252 ed_lam
.grpp
.ener
[egLJSR
],
253 ed_lam
.grpp
.ener
[egCOULSR
],
254 enerd
->grpp
.ener
[egGB
], box_size
,nrnb
,
255 lam_i
,&dvdl_dum
,-1,-1,
256 GMX_DONB_FOREIGNLAMBDA
);
257 sum_epot(&ir
->opts
,&ed_lam
);
258 enerd
->enerpart_lambda
[i
] += ed_lam
.term
[F_EPOT
];
260 destroy_enerdata(&ed_lam
);
264 /* If we are doing GB, calculate bonded forces and apply corrections
265 * to the solvation forces */
266 if (ir
->implicit_solvent
) {
267 calc_gb_forces(cr
,md
,born
,top
,atype
,x
,f
,fr
,idef
,
268 ir
->gb_algorithm
,ir
->sa_algorithm
,nrnb
,bBornRadii
,&pbc
,graph
,enerd
);
279 if (ir
->sc_alpha
!= 0)
281 enerd
->dvdl_nonlin
+= dvdlambda
;
285 enerd
->dvdl_lin
+= dvdlambda
;
290 for(i
=0; i
<enerd
->grpp
.nener
; i
++)
294 enerd
->grpp
.ener
[egBHAMSR
][i
] :
295 enerd
->grpp
.ener
[egLJSR
][i
])
296 + enerd
->grpp
.ener
[egCOULSR
][i
] + enerd
->grpp
.ener
[egGB
][i
];
299 PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",Vsr
,dvdlambda
);
302 GMX_MPE_LOG(ev_do_fnbf_finish
);
306 pr_rvecs(debug
,0,"fshift after SR",fr
->fshift
,SHIFTS
);
309 /* Shift the coordinates. Must be done before bonded forces and PPPM,
310 * but is also necessary for SHAKE and update, therefore it can NOT
311 * go when no bonded forces have to be evaluated.
314 /* Here sometimes we would not need to shift with NBFonly,
315 * but we do so anyhow for consistency of the returned coordinates.
319 shift_self(graph
,box
,x
);
322 inc_nrnb(nrnb
,eNR_SHIFTX
,2*graph
->nnodes
);
326 inc_nrnb(nrnb
,eNR_SHIFTX
,graph
->nnodes
);
329 /* Check whether we need to do bondeds or correct for exclusions */
331 ((flags
& GMX_FORCE_BONDED
)
332 || EEL_RF(fr
->eeltype
) || EEL_FULL(fr
->eeltype
)))
334 /* Since all atoms are in the rectangular or triclinic unit-cell,
335 * only single box vector shifts (2 in x) are required.
337 set_pbc_dd(&pbc
,fr
->ePBC
,cr
->dd
,TRUE
,box
);
341 if (flags
& GMX_FORCE_BONDED
)
343 GMX_MPE_LOG(ev_calc_bonds_start
);
344 calc_bonds(fplog
,cr
->ms
,
345 idef
,x
,hist
,f
,fr
,&pbc
,graph
,enerd
,nrnb
,lambda
,md
,fcd
,
346 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: NULL
, atype
, born
,
347 fr
->bSepDVDL
&& do_per_step(step
,ir
->nstlog
),step
);
349 /* Check if we have to determine energy differences
350 * at foreign lambda's.
352 if (ir
->n_flambda
> 0 && (flags
& GMX_FORCE_DHDL
) &&
353 idef
->ilsort
!= ilsortNO_FE
)
355 if (idef
->ilsort
!= ilsortFE_SORTED
)
357 gmx_incons("The bonded interactions are not sorted for free energy");
359 init_enerdata(mtop
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,&ed_lam
);
361 for(i
=0; i
<enerd
->n_lambda
; i
++)
363 lam_i
= (i
==0 ? lambda
: ir
->flambda
[i
-1]);
365 reset_enerdata(&ir
->opts
,fr
,TRUE
,&ed_lam
,FALSE
);
366 calc_bonds_lambda(fplog
,
367 idef
,x
,fr
,&pbc
,graph
,&ed_lam
,nrnb
,lam_i
,md
,
369 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: NULL
);
370 sum_epot(&ir
->opts
,&ed_lam
);
371 enerd
->enerpart_lambda
[i
] += ed_lam
.term
[F_EPOT
];
373 destroy_enerdata(&ed_lam
);
376 GMX_MPE_LOG(ev_calc_bonds_finish
);
382 if (EEL_FULL(fr
->eeltype
))
384 bSB
= (ir
->nwall
== 2);
388 svmul(ir
->wall_ewald_zfac
,boxs
[ZZ
],boxs
[ZZ
]);
389 box_size
[ZZ
] *= ir
->wall_ewald_zfac
;
392 clear_mat(fr
->vir_el_recip
);
399 Vcorr
= ewald_LRcorrection(fplog
,md
->start
,md
->start
+md
->homenr
,
402 md
->nChargePerturbed
? md
->chargeB
: NULL
,
403 excl
,x
,bSB
? boxs
: box
,mu_tot
,
406 lambda
,&dvdlambda
,&vdip
,&vcharge
);
407 PRINT_SEPDVDL("Ewald excl./charge/dip. corr.",Vcorr
,dvdlambda
);
408 enerd
->dvdl_lin
+= dvdlambda
;
412 if (ir
->ewald_geometry
!= eewg3D
|| ir
->epsilon_surface
!= 0)
414 gmx_fatal(FARGS
,"TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
416 /* The TPI molecule does not have exclusions with the rest
417 * of the system and no intra-molecular PME grid contributions
418 * will be calculated in gmx_pme_calc_energy.
431 case eelPMEUSERSWITCH
:
433 if (cr
->duty
& DUTY_PME
)
435 assert(fr
->n_tpi
>= 0);
436 if (fr
->n_tpi
== 0 || (flags
& GMX_FORCE_STATECHANGED
))
438 pme_flags
= GMX_PME_SPREAD_Q
| GMX_PME_SOLVE
;
439 if (flags
& GMX_FORCE_FORCES
)
441 pme_flags
|= GMX_PME_CALC_F
;
443 if (flags
& GMX_FORCE_VIRIAL
)
445 pme_flags
|= GMX_PME_CALC_ENER_VIR
;
449 /* We don't calculate f, but we do want the potential */
450 pme_flags
|= GMX_PME_CALC_POT
;
452 wallcycle_start(wcycle
,ewcPMEMESH
);
453 status
= gmx_pme_do(fr
->pmedata
,
454 md
->start
,md
->homenr
- fr
->n_tpi
,
456 md
->chargeA
,md
->chargeB
,
458 DOMAINDECOMP(cr
) ? dd_pme_maxshift_x(cr
->dd
) : 0,
459 DOMAINDECOMP(cr
) ? dd_pme_maxshift_y(cr
->dd
) : 0,
461 fr
->vir_el_recip
,fr
->ewaldcoeff
,
462 &Vlr
,lambda
,&dvdlambda
,
464 *cycles_pme
= wallcycle_stop(wcycle
,ewcPMEMESH
);
466 /* We should try to do as little computation after
467 * this as possible, because parallel PME synchronizes
468 * the nodes, so we want all load imbalance of the rest
469 * of the force calculation to be before the PME call.
470 * DD load balancing is done on the whole time of
471 * the force call (without PME).
476 /* Determine the PME grid energy of the test molecule
477 * with the PME grid potential of the other charges.
479 gmx_pme_calc_energy(fr
->pmedata
,fr
->n_tpi
,
480 x
+ md
->homenr
- fr
->n_tpi
,
481 md
->chargeA
+ md
->homenr
- fr
->n_tpi
,
484 PRINT_SEPDVDL("PME mesh",Vlr
,dvdlambda
);
488 /* Energies and virial are obtained later from the PME nodes */
489 /* but values have to be zeroed out here */
494 Vlr
= do_ewald(fplog
,FALSE
,ir
,x
,fr
->f_novirsum
,
495 md
->chargeA
,md
->chargeB
,
496 box_size
,cr
,md
->homenr
,
497 fr
->vir_el_recip
,fr
->ewaldcoeff
,
498 lambda
,&dvdlambda
,fr
->ewald_table
);
499 PRINT_SEPDVDL("Ewald long-range",Vlr
,dvdlambda
);
503 gmx_fatal(FARGS
,"No such electrostatics method implemented %s",
504 eel_names
[fr
->eeltype
]);
508 gmx_fatal(FARGS
,"Error %d in long range electrostatics routine %s",
509 status
,EELTYPE(fr
->eeltype
));
511 enerd
->dvdl_lin
+= dvdlambda
;
512 enerd
->term
[F_COUL_RECIP
] = Vlr
+ Vcorr
;
515 fprintf(debug
,"Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
516 Vlr
,Vcorr
,enerd
->term
[F_COUL_RECIP
]);
517 pr_rvecs(debug
,0,"vir_el_recip after corr",fr
->vir_el_recip
,DIM
);
518 pr_rvecs(debug
,0,"fshift after LR Corrections",fr
->fshift
,SHIFTS
);
523 if (EEL_RF(fr
->eeltype
))
527 if (fr
->eeltype
!= eelRF_NEC
)
529 enerd
->term
[F_RF_EXCL
] =
530 RF_excl_correction(fplog
,fr
,graph
,md
,excl
,x
,f
,
531 fr
->fshift
,&pbc
,lambda
,&dvdlambda
);
534 enerd
->dvdl_lin
+= dvdlambda
;
535 PRINT_SEPDVDL("RF exclusion correction",
536 enerd
->term
[F_RF_EXCL
],dvdlambda
);
544 print_nrnb(debug
,nrnb
);
552 MPI_Barrier(cr
->mpi_comm_mygroup
);
555 if (fr
->timesteps
== 11)
557 fprintf(stderr
,"* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
558 cr
->nodeid
, gmx_step_str(fr
->timesteps
,buf
),
559 100*fr
->t_wait
/(fr
->t_wait
+fr
->t_fnbf
),
560 (fr
->t_fnbf
+fr
->t_wait
)/fr
->t_fnbf
);
568 pr_rvecs(debug
,0,"fshift after bondeds",fr
->fshift
,SHIFTS
);
571 GMX_MPE_LOG(ev_force_finish
);
575 void init_enerdata(int ngener
,int n_flambda
,gmx_enerdata_t
*enerd
)
579 for(i
=0; i
<F_NRE
; i
++)
587 fprintf(debug
,"Creating %d sized group matrix for energies\n",n2
);
589 enerd
->grpp
.nener
= n2
;
590 for(i
=0; (i
<egNR
); i
++)
592 snew(enerd
->grpp
.ener
[i
],n2
);
597 enerd
->n_lambda
= 1 + n_flambda
;
598 snew(enerd
->enerpart_lambda
,enerd
->n_lambda
);
606 void destroy_enerdata(gmx_enerdata_t
*enerd
)
610 for(i
=0; (i
<egNR
); i
++)
612 sfree(enerd
->grpp
.ener
[i
]);
617 sfree(enerd
->enerpart_lambda
);
621 static real
sum_v(int n
,real v
[])
633 void sum_epot(t_grpopts
*opts
,gmx_enerdata_t
*enerd
)
635 gmx_grppairener_t
*grpp
;
642 /* Accumulate energies */
643 epot
[F_COUL_SR
] = sum_v(grpp
->nener
,grpp
->ener
[egCOULSR
]);
644 epot
[F_LJ
] = sum_v(grpp
->nener
,grpp
->ener
[egLJSR
]);
645 epot
[F_LJ14
] = sum_v(grpp
->nener
,grpp
->ener
[egLJ14
]);
646 epot
[F_COUL14
] = sum_v(grpp
->nener
,grpp
->ener
[egCOUL14
]);
647 epot
[F_COUL_LR
] = sum_v(grpp
->nener
,grpp
->ener
[egCOULLR
]);
648 epot
[F_LJ_LR
] = sum_v(grpp
->nener
,grpp
->ener
[egLJLR
]);
649 /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
650 epot
[F_GBPOL
] += sum_v(grpp
->nener
,grpp
->ener
[egGB
]);
652 /* lattice part of LR doesnt belong to any group
653 * and has been added earlier
655 epot
[F_BHAM
] = sum_v(grpp
->nener
,grpp
->ener
[egBHAMSR
]);
656 epot
[F_BHAM_LR
] = sum_v(grpp
->nener
,grpp
->ener
[egBHAMLR
]);
659 for(i
=0; (i
<F_EPOT
); i
++)
660 if (i
!= F_DISRESVIOL
&& i
!= F_ORIRESDEV
&& i
!= F_DIHRESVIOL
)
661 epot
[F_EPOT
] += epot
[i
];
664 void sum_dhdl(gmx_enerdata_t
*enerd
,double lambda
,t_inputrec
*ir
)
667 double dlam
,dhdl_lin
;
669 enerd
->term
[F_DVDL
] = enerd
->dvdl_lin
+ enerd
->dvdl_nonlin
;
673 fprintf(debug
,"dvdl: %f, non-linear %f + linear %f\n",
674 enerd
->term
[F_DVDL
],enerd
->dvdl_nonlin
,enerd
->dvdl_lin
);
677 /* Notes on the foreign lambda free energy difference evaluation:
678 * Adding the potential and ekin terms that depend linearly on lambda
679 * as delta lam * dvdl to the energy differences is exact.
680 * For the constraint dvdl this is not exact, but we have no other option.
681 * For the non-bonded LR term we assume that the soft-core (if present)
682 * no longer affects the energy beyond the short-range cut-off,
683 * which is a very good approximation (except for exotic settings).
685 for(i
=1; i
<enerd
->n_lambda
; i
++)
687 dlam
= (ir
->flambda
[i
-1] - lambda
);
689 enerd
->dvdl_lin
+ enerd
->term
[F_DKDL
] + enerd
->term
[F_DHDL_CON
];
692 fprintf(debug
,"enerdiff lam %g: non-linear %f linear %f*%f\n",
694 enerd
->enerpart_lambda
[i
] - enerd
->enerpart_lambda
[0],
697 enerd
->enerpart_lambda
[i
] += dlam
*dhdl_lin
;
702 void reset_enerdata(t_grpopts
*opts
,
703 t_forcerec
*fr
,gmx_bool bNS
,
704 gmx_enerdata_t
*enerd
,
710 /* First reset all energy components, except for the long range terms
711 * on the master at non neighbor search steps, since the long range
712 * terms have already been summed at the last neighbor search step.
714 bKeepLR
= (fr
->bTwinRange
&& !bNS
);
715 for(i
=0; (i
<egNR
); i
++) {
716 if (!(bKeepLR
&& bMaster
&& (i
== egCOULLR
|| i
== egLJLR
))) {
717 for(j
=0; (j
<enerd
->grpp
.nener
); j
++)
718 enerd
->grpp
.ener
[i
][j
] = 0.0;
721 enerd
->dvdl_lin
= 0.0;
722 enerd
->dvdl_nonlin
= 0.0;
724 /* Normal potential energy components */
725 for(i
=0; (i
<=F_EPOT
); i
++) {
726 enerd
->term
[i
] = 0.0;
728 /* Initialize the dVdlambda term with the long range contribution */
729 enerd
->term
[F_DVDL
] = 0.0;
730 enerd
->term
[F_DKDL
] = 0.0;
731 enerd
->term
[F_DHDL_CON
] = 0.0;
732 if (enerd
->n_lambda
> 0)
734 for(i
=0; i
<enerd
->n_lambda
; i
++)
736 enerd
->enerpart_lambda
[i
] = 0.0;