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
49 #include "nonbonded.h"
63 #include "mtop_util.h"
67 /* MSVC definition for __cpuid() */
73 t_forcerec
*mk_forcerec(void)
83 static void pr_nbfp(FILE *fp
,real
*nbfp
,gmx_bool bBHAM
,int atnr
)
87 for(i
=0; (i
<atnr
); i
++) {
88 for(j
=0; (j
<atnr
); j
++) {
89 fprintf(fp
,"%2d - %2d",i
,j
);
91 fprintf(fp
," a=%10g, b=%10g, c=%10g\n",BHAMA(nbfp
,atnr
,i
,j
),
92 BHAMB(nbfp
,atnr
,i
,j
),BHAMC(nbfp
,atnr
,i
,j
));
94 fprintf(fp
," c6=%10g, c12=%10g\n",C6(nbfp
,atnr
,i
,j
),
101 static real
*mk_nbfp(const gmx_ffparams_t
*idef
,gmx_bool bBHAM
)
108 snew(nbfp
,3*atnr
*atnr
);
109 for(i
=k
=0; (i
<atnr
); i
++) {
110 for(j
=0; (j
<atnr
); j
++,k
++) {
111 BHAMA(nbfp
,atnr
,i
,j
) = idef
->iparams
[k
].bham
.a
;
112 BHAMB(nbfp
,atnr
,i
,j
) = idef
->iparams
[k
].bham
.b
;
113 BHAMC(nbfp
,atnr
,i
,j
) = idef
->iparams
[k
].bham
.c
;
118 snew(nbfp
,2*atnr
*atnr
);
119 for(i
=k
=0; (i
<atnr
); i
++) {
120 for(j
=0; (j
<atnr
); j
++,k
++) {
121 C6(nbfp
,atnr
,i
,j
) = idef
->iparams
[k
].lj
.c6
;
122 C12(nbfp
,atnr
,i
,j
) = idef
->iparams
[k
].lj
.c12
;
129 /* This routine sets fr->solvent_opt to the most common solvent in the
130 * system, e.g. esolSPC or esolTIP4P. It will also mark each charge group in
131 * the fr->solvent_type array with the correct type (or esolNO).
133 * Charge groups that fulfill the conditions but are not identical to the
134 * most common one will be marked as esolNO in the solvent_type array.
136 * TIP3p is identical to SPC for these purposes, so we call it
137 * SPC in the arrays (Apologies to Bill Jorgensen ;-)
139 * NOTE: QM particle should not
140 * become an optimized solvent. Not even if there is only one charge
150 } solvent_parameters_t
;
153 check_solvent_cg(const gmx_moltype_t
*molt
,
156 const unsigned char *qm_grpnr
,
157 const t_grps
*qm_grps
,
159 int *n_solvent_parameters
,
160 solvent_parameters_t
**solvent_parameters_p
,
164 const t_blocka
* excl
;
175 solvent_parameters_t
*solvent_parameters
;
177 /* We use a list with parameters for each solvent type.
178 * Every time we discover a new molecule that fulfills the basic
179 * conditions for a solvent we compare with the previous entries
180 * in these lists. If the parameters are the same we just increment
181 * the counter for that type, and otherwise we create a new type
182 * based on the current molecule.
184 * Once we've finished going through all molecules we check which
185 * solvent is most common, and mark all those molecules while we
186 * clear the flag on all others.
189 solvent_parameters
= *solvent_parameters_p
;
191 /* Mark the cg first as non optimized */
194 /* Check if this cg has no exclusions with atoms in other charge groups
195 * and all atoms inside the charge group excluded.
196 * We only have 3 or 4 atom solvent loops.
198 if (GET_CGINFO_EXCL_INTER(cginfo
) ||
199 !GET_CGINFO_EXCL_INTRA(cginfo
))
204 /* Get the indices of the first atom in this charge group */
205 j0
= molt
->cgs
.index
[cg0
];
206 j1
= molt
->cgs
.index
[cg0
+1];
208 /* Number of atoms in our molecule */
213 "Moltype '%s': there are %d atoms in this charge group\n",
217 /* Check if it could be an SPC (3 atoms) or TIP4p (4) water,
225 /* Check if we are doing QM on this group */
227 if (qm_grpnr
!= NULL
)
229 for(j
=j0
; j
<j1
&& !qm
; j
++)
231 qm
= (qm_grpnr
[j
] < qm_grps
->nr
- 1);
234 /* Cannot use solvent optimization with QM */
240 atom
= molt
->atoms
.atom
;
242 /* Still looks like a solvent, time to check parameters */
244 /* If it is perturbed (free energy) we can't use the solvent loops,
245 * so then we just skip to the next molecule.
249 for(j
=j0
; j
<j1
&& !perturbed
; j
++)
251 perturbed
= PERTURBED(atom
[j
]);
259 /* Now it's only a question if the VdW and charge parameters
260 * are OK. Before doing the check we compare and see if they are
261 * identical to a possible previous solvent type.
262 * First we assign the current types and charges.
266 tmp_vdwtype
[j
] = atom
[j0
+j
].type
;
267 tmp_charge
[j
] = atom
[j0
+j
].q
;
270 /* Does it match any previous solvent type? */
271 for(k
=0 ; k
<*n_solvent_parameters
; k
++)
276 /* We can only match SPC with 3 atoms and TIP4p with 4 atoms */
277 if( (solvent_parameters
[k
].model
==esolSPC
&& nj
!=3) ||
278 (solvent_parameters
[k
].model
==esolTIP4P
&& nj
!=4) )
281 /* Check that types & charges match for all atoms in molecule */
282 for(j
=0 ; j
<nj
&& match
==TRUE
; j
++)
284 if (tmp_vdwtype
[j
] != solvent_parameters
[k
].vdwtype
[j
])
288 if(tmp_charge
[j
] != solvent_parameters
[k
].charge
[j
])
295 /* Congratulations! We have a matched solvent.
296 * Flag it with this type for later processing.
299 solvent_parameters
[k
].count
+= nmol
;
301 /* We are done with this charge group */
306 /* If we get here, we have a tentative new solvent type.
307 * Before we add it we must check that it fulfills the requirements
308 * of the solvent optimized loops. First determine which atoms have
314 tjA
= tmp_vdwtype
[j
];
316 /* Go through all other tpes and see if any have non-zero
317 * VdW parameters when combined with this one.
319 for(k
=0; k
<fr
->ntype
&& (has_vdw
[j
]==FALSE
); k
++)
321 /* We already checked that the atoms weren't perturbed,
322 * so we only need to check state A now.
326 has_vdw
[j
] = (has_vdw
[j
] ||
327 (BHAMA(fr
->nbfp
,fr
->ntype
,tjA
,k
) != 0.0) ||
328 (BHAMB(fr
->nbfp
,fr
->ntype
,tjA
,k
) != 0.0) ||
329 (BHAMC(fr
->nbfp
,fr
->ntype
,tjA
,k
) != 0.0));
334 has_vdw
[j
] = (has_vdw
[j
] ||
335 (C6(fr
->nbfp
,fr
->ntype
,tjA
,k
) != 0.0) ||
336 (C12(fr
->nbfp
,fr
->ntype
,tjA
,k
) != 0.0));
341 /* Now we know all we need to make the final check and assignment. */
345 * For this we require thatn all atoms have charge,
346 * the charges on atom 2 & 3 should be the same, and only
347 * atom 1 should have VdW.
349 if (has_vdw
[0] == TRUE
&&
350 has_vdw
[1] == FALSE
&&
351 has_vdw
[2] == FALSE
&&
352 tmp_charge
[0] != 0 &&
353 tmp_charge
[1] != 0 &&
354 tmp_charge
[2] == tmp_charge
[1])
356 srenew(solvent_parameters
,*n_solvent_parameters
+1);
357 solvent_parameters
[*n_solvent_parameters
].model
= esolSPC
;
358 solvent_parameters
[*n_solvent_parameters
].count
= nmol
;
361 solvent_parameters
[*n_solvent_parameters
].vdwtype
[k
] = tmp_vdwtype
[k
];
362 solvent_parameters
[*n_solvent_parameters
].charge
[k
] = tmp_charge
[k
];
365 *cg_sp
= *n_solvent_parameters
;
366 (*n_solvent_parameters
)++;
371 /* Or could it be a TIP4P?
372 * For this we require thatn atoms 2,3,4 have charge, but not atom 1.
373 * Only atom 1 should have VdW.
375 if(has_vdw
[0] == TRUE
&&
376 has_vdw
[1] == FALSE
&&
377 has_vdw
[2] == FALSE
&&
378 has_vdw
[3] == FALSE
&&
379 tmp_charge
[0] == 0 &&
380 tmp_charge
[1] != 0 &&
381 tmp_charge
[2] == tmp_charge
[1] &&
384 srenew(solvent_parameters
,*n_solvent_parameters
+1);
385 solvent_parameters
[*n_solvent_parameters
].model
= esolTIP4P
;
386 solvent_parameters
[*n_solvent_parameters
].count
= nmol
;
389 solvent_parameters
[*n_solvent_parameters
].vdwtype
[k
] = tmp_vdwtype
[k
];
390 solvent_parameters
[*n_solvent_parameters
].charge
[k
] = tmp_charge
[k
];
393 *cg_sp
= *n_solvent_parameters
;
394 (*n_solvent_parameters
)++;
398 *solvent_parameters_p
= solvent_parameters
;
402 check_solvent(FILE * fp
,
403 const gmx_mtop_t
* mtop
,
405 cginfo_mb_t
*cginfo_mb
)
408 const t_block
* mols
;
409 const gmx_moltype_t
*molt
;
410 int mb
,mol
,cg_mol
,at_offset
,cg_offset
,am
,cgm
,i
,nmol_ch
,nmol
;
411 int n_solvent_parameters
;
412 solvent_parameters_t
*solvent_parameters
;
418 fprintf(debug
,"Going to determine what solvent types we have.\n");
423 n_solvent_parameters
= 0;
424 solvent_parameters
= NULL
;
425 /* Allocate temporary array for solvent type */
426 snew(cg_sp
,mtop
->nmolblock
);
430 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
432 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
434 /* Here we have to loop over all individual molecules
435 * because we need to check for QMMM particles.
437 snew(cg_sp
[mb
],cginfo_mb
[mb
].cg_mod
);
438 nmol_ch
= cginfo_mb
[mb
].cg_mod
/cgs
->nr
;
439 nmol
= mtop
->molblock
[mb
].nmol
/nmol_ch
;
440 for(mol
=0; mol
<nmol_ch
; mol
++)
443 am
= mol
*cgs
->index
[cgs
->nr
];
444 for(cg_mol
=0; cg_mol
<cgs
->nr
; cg_mol
++)
446 check_solvent_cg(molt
,cg_mol
,nmol
,
447 mtop
->groups
.grpnr
[egcQMMM
] ?
448 mtop
->groups
.grpnr
[egcQMMM
]+at_offset
+am
: 0,
449 &mtop
->groups
.grps
[egcQMMM
],
451 &n_solvent_parameters
,&solvent_parameters
,
452 cginfo_mb
[mb
].cginfo
[cgm
+cg_mol
],
453 &cg_sp
[mb
][cgm
+cg_mol
]);
456 cg_offset
+= cgs
->nr
;
457 at_offset
+= cgs
->index
[cgs
->nr
];
460 /* Puh! We finished going through all charge groups.
461 * Now find the most common solvent model.
464 /* Most common solvent this far */
466 for(i
=0;i
<n_solvent_parameters
;i
++)
469 solvent_parameters
[i
].count
> solvent_parameters
[bestsp
].count
)
477 bestsol
= solvent_parameters
[bestsp
].model
;
484 #ifdef DISABLE_WATER_NLIST
489 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
491 cgs
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].cgs
;
492 nmol
= (mtop
->molblock
[mb
].nmol
*cgs
->nr
)/cginfo_mb
[mb
].cg_mod
;
493 for(i
=0; i
<cginfo_mb
[mb
].cg_mod
; i
++)
495 if (cg_sp
[mb
][i
] == bestsp
)
497 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[i
],bestsol
);
502 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[i
],esolNO
);
509 if (bestsol
!= esolNO
&& fp
!=NULL
)
511 fprintf(fp
,"\nEnabling %s-like water optimization for %d molecules.\n\n",
513 solvent_parameters
[bestsp
].count
);
516 sfree(solvent_parameters
);
517 fr
->solvent_opt
= bestsol
;
520 static cginfo_mb_t
*init_cginfo_mb(FILE *fplog
,const gmx_mtop_t
*mtop
,
521 t_forcerec
*fr
,gmx_bool bNoSolvOpt
)
524 const t_blocka
*excl
;
525 const gmx_moltype_t
*molt
;
526 const gmx_molblock_t
*molb
;
527 cginfo_mb_t
*cginfo_mb
;
529 int cg_offset
,a_offset
,cgm
,am
;
530 int mb
,m
,ncg_tot
,cg
,a0
,a1
,gid
,ai
,j
,aj
,excl_nalloc
;
531 gmx_bool bId
,*bExcl
,bExclIntraAll
,bExclInter
;
533 ncg_tot
= ncg_mtop(mtop
);
534 snew(cginfo_mb
,mtop
->nmolblock
);
537 snew(bExcl
,excl_nalloc
);
540 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
542 molb
= &mtop
->molblock
[mb
];
543 molt
= &mtop
->moltype
[molb
->type
];
547 /* Check if the cginfo is identical for all molecules in this block.
548 * If so, we only need an array of the size of one molecule.
549 * Otherwise we make an array of #mol times #cgs per molecule.
553 for(m
=0; m
<molb
->nmol
; m
++)
555 am
= m
*cgs
->index
[cgs
->nr
];
556 for(cg
=0; cg
<cgs
->nr
; cg
++)
559 a1
= cgs
->index
[cg
+1];
560 if (ggrpnr(&mtop
->groups
,egcENER
,a_offset
+am
+a0
) !=
561 ggrpnr(&mtop
->groups
,egcENER
,a_offset
+a0
))
565 if (mtop
->groups
.grpnr
[egcQMMM
] != NULL
)
567 for(ai
=a0
; ai
<a1
; ai
++)
569 if (mtop
->groups
.grpnr
[egcQMMM
][a_offset
+am
+ai
] !=
570 mtop
->groups
.grpnr
[egcQMMM
][a_offset
+ai
])
579 cginfo_mb
[mb
].cg_start
= cg_offset
;
580 cginfo_mb
[mb
].cg_end
= cg_offset
+ molb
->nmol
*cgs
->nr
;
581 cginfo_mb
[mb
].cg_mod
= (bId
? 1 : molb
->nmol
)*cgs
->nr
;
582 snew(cginfo_mb
[mb
].cginfo
,cginfo_mb
[mb
].cg_mod
);
583 cginfo
= cginfo_mb
[mb
].cginfo
;
585 for(m
=0; m
<(bId
? 1 : molb
->nmol
); m
++)
588 am
= m
*cgs
->index
[cgs
->nr
];
589 for(cg
=0; cg
<cgs
->nr
; cg
++)
592 a1
= cgs
->index
[cg
+1];
594 /* Store the energy group in cginfo */
595 gid
= ggrpnr(&mtop
->groups
,egcENER
,a_offset
+am
+a0
);
596 SET_CGINFO_GID(cginfo
[cgm
+cg
],gid
);
598 /* Check the intra/inter charge group exclusions */
599 if (a1
-a0
> excl_nalloc
) {
600 excl_nalloc
= a1
- a0
;
601 srenew(bExcl
,excl_nalloc
);
603 /* bExclIntraAll: all intra cg interactions excluded
604 * bExclInter: any inter cg interactions excluded
606 bExclIntraAll
= TRUE
;
608 for(ai
=a0
; ai
<a1
; ai
++) {
609 /* Clear the exclusion list for atom ai */
610 for(aj
=a0
; aj
<a1
; aj
++) {
611 bExcl
[aj
-a0
] = FALSE
;
613 /* Loop over all the exclusions of atom ai */
614 for(j
=excl
->index
[ai
]; j
<excl
->index
[ai
+1]; j
++)
617 if (aj
< a0
|| aj
>= a1
)
626 /* Check if ai excludes a0 to a1 */
627 for(aj
=a0
; aj
<a1
; aj
++)
631 bExclIntraAll
= FALSE
;
637 SET_CGINFO_EXCL_INTRA(cginfo
[cgm
+cg
]);
641 SET_CGINFO_EXCL_INTER(cginfo
[cgm
+cg
]);
643 if (a1
- a0
> MAX_CHARGEGROUP_SIZE
)
645 /* The size in cginfo is currently only read with DD */
646 gmx_fatal(FARGS
,"A charge group has size %d which is larger than the limit of %d atoms",a1
-a0
,MAX_CHARGEGROUP_SIZE
);
648 SET_CGINFO_NATOMS(cginfo
[cgm
+cg
],a1
-a0
);
651 cg_offset
+= molb
->nmol
*cgs
->nr
;
652 a_offset
+= molb
->nmol
*cgs
->index
[cgs
->nr
];
656 /* the solvent optimizer is called after the QM is initialized,
657 * because we don't want to have the QM subsystemto become an
661 check_solvent(fplog
,mtop
,fr
,cginfo_mb
);
663 if (getenv("GMX_NO_SOLV_OPT"))
667 fprintf(fplog
,"Found environment variable GMX_NO_SOLV_OPT.\n"
668 "Disabling all solvent optimization\n");
670 fr
->solvent_opt
= esolNO
;
674 fr
->solvent_opt
= esolNO
;
676 if (!fr
->solvent_opt
)
678 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
680 for(cg
=0; cg
<cginfo_mb
[mb
].cg_mod
; cg
++)
682 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[cg
],esolNO
);
690 static int *cginfo_expand(int nmb
,cginfo_mb_t
*cgi_mb
)
695 ncg
= cgi_mb
[nmb
-1].cg_end
;
698 for(cg
=0; cg
<ncg
; cg
++)
700 while (cg
>= cgi_mb
[mb
].cg_end
)
705 cgi_mb
[mb
].cginfo
[(cg
- cgi_mb
[mb
].cg_start
) % cgi_mb
[mb
].cg_mod
];
711 static void set_chargesum(FILE *log
,t_forcerec
*fr
,const gmx_mtop_t
*mtop
)
715 const t_atoms
*atoms
;
718 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
720 nmol
= mtop
->molblock
[mb
].nmol
;
721 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
722 for(i
=0; i
<atoms
->nr
; i
++)
724 qsum
+= nmol
*atoms
->atom
[i
].q
;
728 if (fr
->efep
!= efepNO
)
731 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
733 nmol
= mtop
->molblock
[mb
].nmol
;
734 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
735 for(i
=0; i
<atoms
->nr
; i
++)
737 qsum
+= nmol
*atoms
->atom
[i
].qB
;
744 fr
->qsum
[1] = fr
->qsum
[0];
747 if (fr
->efep
== efepNO
)
748 fprintf(log
,"System total charge: %.3f\n",fr
->qsum
[0]);
750 fprintf(log
,"System total charge, top. A: %.3f top. B: %.3f\n",
751 fr
->qsum
[0],fr
->qsum
[1]);
755 void update_forcerec(FILE *log
,t_forcerec
*fr
,matrix box
)
757 if (fr
->eeltype
== eelGRF
)
759 calc_rffac(NULL
,fr
->eeltype
,fr
->epsilon_r
,fr
->epsilon_rf
,
760 fr
->rcoulomb
,fr
->temp
,fr
->zsquare
,box
,
761 &fr
->kappa
,&fr
->k_rf
,&fr
->c_rf
);
765 void set_avcsixtwelve(FILE *fplog
,t_forcerec
*fr
,const gmx_mtop_t
*mtop
)
767 const t_atoms
*atoms
,*atoms_tpi
;
768 const t_blocka
*excl
;
769 int mb
,nmol
,nmolc
,i
,j
,tpi
,tpj
,j1
,j2
,k
,n
,nexcl
,q
;
770 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
771 long long int npair
,npair_ij
,tmpi
,tmpj
;
773 double npair
, npair_ij
,tmpi
,tmpj
;
784 for(q
=0; q
<(fr
->efep
==efepNO
? 1 : 2); q
++) {
790 /* Count the types so we avoid natoms^2 operations */
792 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
793 nmol
= mtop
->molblock
[mb
].nmol
;
794 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
795 for(i
=0; i
<atoms
->nr
; i
++) {
798 tpi
= atoms
->atom
[i
].type
;
802 tpi
= atoms
->atom
[i
].typeB
;
804 typecount
[tpi
] += nmol
;
807 for(tpi
=0; tpi
<ntp
; tpi
++) {
808 for(tpj
=tpi
; tpj
<ntp
; tpj
++) {
809 tmpi
= typecount
[tpi
];
810 tmpj
= typecount
[tpj
];
813 npair_ij
= tmpi
*tmpj
;
817 npair_ij
= tmpi
*(tmpi
- 1)/2;
820 csix
+= npair_ij
*BHAMC(nbfp
,ntp
,tpi
,tpj
);
822 csix
+= npair_ij
* C6(nbfp
,ntp
,tpi
,tpj
);
823 ctwelve
+= npair_ij
* C12(nbfp
,ntp
,tpi
,tpj
);
829 /* Subtract the excluded pairs.
830 * The main reason for substracting exclusions is that in some cases
831 * some combinations might never occur and the parameters could have
832 * any value. These unused values should not influence the dispersion
835 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
836 nmol
= mtop
->molblock
[mb
].nmol
;
837 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
838 excl
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].excls
;
839 for(i
=0; (i
<atoms
->nr
); i
++) {
842 tpi
= atoms
->atom
[i
].type
;
846 tpi
= atoms
->atom
[i
].typeB
;
849 j2
= excl
->index
[i
+1];
850 for(j
=j1
; j
<j2
; j
++) {
856 tpj
= atoms
->atom
[k
].type
;
860 tpj
= atoms
->atom
[k
].typeB
;
863 csix
-= nmol
*BHAMC(nbfp
,ntp
,tpi
,tpj
);
865 csix
-= nmol
*C6 (nbfp
,ntp
,tpi
,tpj
);
866 ctwelve
-= nmol
*C12(nbfp
,ntp
,tpi
,tpj
);
874 /* Only correct for the interaction of the test particle
875 * with the rest of the system.
878 &mtop
->moltype
[mtop
->molblock
[mtop
->nmolblock
-1].type
].atoms
;
881 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
882 nmol
= mtop
->molblock
[mb
].nmol
;
883 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
884 for(j
=0; j
<atoms
->nr
; j
++) {
886 /* Remove the interaction of the test charge group
889 if (mb
== mtop
->nmolblock
-1)
893 if (mb
== 0 && nmol
== 1)
895 gmx_fatal(FARGS
,"Old format tpr with TPI, please generate a new tpr file");
900 tpj
= atoms
->atom
[j
].type
;
904 tpj
= atoms
->atom
[j
].typeB
;
906 for(i
=0; i
<fr
->n_tpi
; i
++)
910 tpi
= atoms_tpi
->atom
[i
].type
;
914 tpi
= atoms_tpi
->atom
[i
].typeB
;
918 csix
+= nmolc
*BHAMC(nbfp
,ntp
,tpi
,tpj
);
922 csix
+= nmolc
*C6 (nbfp
,ntp
,tpi
,tpj
);
923 ctwelve
+= nmolc
*C12(nbfp
,ntp
,tpi
,tpj
);
930 if (npair
- nexcl
<= 0 && fplog
) {
931 fprintf(fplog
,"\nWARNING: There are no atom pairs for dispersion correction\n\n");
935 csix
/= npair
- nexcl
;
936 ctwelve
/= npair
- nexcl
;
939 fprintf(debug
,"Counted %d exclusions\n",nexcl
);
940 fprintf(debug
,"Average C6 parameter is: %10g\n",(double)csix
);
941 fprintf(debug
,"Average C12 parameter is: %10g\n",(double)ctwelve
);
943 fr
->avcsix
[q
] = csix
;
944 fr
->avctwelve
[q
] = ctwelve
;
948 if (fr
->eDispCorr
== edispcAllEner
||
949 fr
->eDispCorr
== edispcAllEnerPres
)
951 fprintf(fplog
,"Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
952 fr
->avcsix
[0],fr
->avctwelve
[0]);
956 fprintf(fplog
,"Long Range LJ corr.: <C6> %10.4e\n",fr
->avcsix
[0]);
962 static void set_bham_b_max(FILE *fplog
,t_forcerec
*fr
,
963 const gmx_mtop_t
*mtop
)
965 const t_atoms
*at1
,*at2
;
966 int mt1
,mt2
,i
,j
,tpi
,tpj
,ntypes
;
972 fprintf(fplog
,"Determining largest Buckingham b parameter for table\n");
979 for(mt1
=0; mt1
<mtop
->nmoltype
; mt1
++)
981 at1
= &mtop
->moltype
[mt1
].atoms
;
982 for(i
=0; (i
<at1
->nr
); i
++)
984 tpi
= at1
->atom
[i
].type
;
986 gmx_fatal(FARGS
,"Atomtype[%d] = %d, maximum = %d",i
,tpi
,ntypes
);
988 for(mt2
=mt1
; mt2
<mtop
->nmoltype
; mt2
++)
990 at2
= &mtop
->moltype
[mt2
].atoms
;
991 for(j
=0; (j
<at2
->nr
); j
++) {
992 tpj
= at2
->atom
[j
].type
;
995 gmx_fatal(FARGS
,"Atomtype[%d] = %d, maximum = %d",j
,tpj
,ntypes
);
997 b
= BHAMB(nbfp
,ntypes
,tpi
,tpj
);
998 if (b
> fr
->bham_b_max
)
1002 if ((b
< bmin
) || (bmin
==-1))
1012 fprintf(fplog
,"Buckingham b parameters, min: %g, max: %g\n",
1013 bmin
,fr
->bham_b_max
);
1017 static void make_nbf_tables(FILE *fp
,const output_env_t oenv
,
1018 t_forcerec
*fr
,real rtab
,
1019 const t_commrec
*cr
,
1020 const char *tabfn
,char *eg1
,char *eg2
,
1026 if (tabfn
== NULL
) {
1028 fprintf(debug
,"No table file name passed, can not read table, can not do non-bonded interactions\n");
1032 sprintf(buf
,"%s",tabfn
);
1034 /* Append the two energy group names */
1035 sprintf(buf
+ strlen(tabfn
) - strlen(ftp2ext(efXVG
)) - 1,"_%s_%s.%s",
1036 eg1
,eg2
,ftp2ext(efXVG
));
1037 nbl
->tab
= make_tables(fp
,oenv
,fr
,MASTER(cr
),buf
,rtab
,0);
1038 /* Copy the contents of the table to separate coulomb and LJ tables too,
1039 * to improve cache performance.
1042 /* For performance reasons we want
1043 * the table data to be aligned to 16-byte. The pointer could be freed
1044 * but currently isn't.
1046 snew_aligned(nbl
->vdwtab
,8*(nbl
->tab
.n
+1),16);
1047 snew_aligned(nbl
->coultab
,4*(nbl
->tab
.n
+1),16);
1049 for(i
=0; i
<=nbl
->tab
.n
; i
++) {
1051 nbl
->coultab
[4*i
+j
] = nbl
->tab
.tab
[12*i
+j
];
1053 nbl
->vdwtab
[8*i
+j
] = nbl
->tab
.tab
[12*i
+4+j
];
1057 static void count_tables(int ftype1
,int ftype2
,const gmx_mtop_t
*mtop
,
1058 int *ncount
,int **count
)
1060 const gmx_moltype_t
*molt
;
1062 int mt
,ftype
,stride
,i
,j
,tabnr
;
1064 for(mt
=0; mt
<mtop
->nmoltype
; mt
++)
1066 molt
= &mtop
->moltype
[mt
];
1067 for(ftype
=0; ftype
<F_NRE
; ftype
++)
1069 if (ftype
== ftype1
|| ftype
== ftype2
) {
1070 il
= &molt
->ilist
[ftype
];
1071 stride
= 1 + NRAL(ftype
);
1072 for(i
=0; i
<il
->nr
; i
+=stride
) {
1073 tabnr
= mtop
->ffparams
.iparams
[il
->iatoms
[i
]].tab
.table
;
1075 gmx_fatal(FARGS
,"A bonded table number is smaller than 0: %d\n",tabnr
);
1076 if (tabnr
>= *ncount
) {
1077 srenew(*count
,tabnr
+1);
1078 for(j
=*ncount
; j
<tabnr
+1; j
++)
1089 static bondedtable_t
*make_bonded_tables(FILE *fplog
,
1090 int ftype1
,int ftype2
,
1091 const gmx_mtop_t
*mtop
,
1092 const char *basefn
,const char *tabext
)
1094 int i
,ncount
,*count
;
1102 count_tables(ftype1
,ftype2
,mtop
,&ncount
,&count
);
1106 for(i
=0; i
<ncount
; i
++) {
1108 sprintf(tabfn
,"%s",basefn
);
1109 sprintf(tabfn
+ strlen(basefn
) - strlen(ftp2ext(efXVG
)) - 1,"_%s%d.%s",
1110 tabext
,i
,ftp2ext(efXVG
));
1111 tab
[i
] = make_bonded_table(fplog
,tabfn
,NRAL(ftype1
)-2);
1120 void forcerec_set_ranges(t_forcerec
*fr
,
1121 int ncg_home
,int ncg_force
,
1123 int natoms_force_constr
,int natoms_f_novirsum
)
1128 /* fr->ncg_force is unused in the standard code,
1129 * but it can be useful for modified code dealing with charge groups.
1131 fr
->ncg_force
= ncg_force
;
1132 fr
->natoms_force
= natoms_force
;
1133 fr
->natoms_force_constr
= natoms_force_constr
;
1135 if (fr
->natoms_force_constr
> fr
->nalloc_force
)
1137 fr
->nalloc_force
= over_alloc_dd(fr
->natoms_force_constr
);
1141 srenew(fr
->f_twin
,fr
->nalloc_force
);
1145 if (fr
->bF_NoVirSum
)
1147 fr
->f_novirsum_n
= natoms_f_novirsum
;
1148 if (fr
->f_novirsum_n
> fr
->f_novirsum_nalloc
)
1150 fr
->f_novirsum_nalloc
= over_alloc_dd(fr
->f_novirsum_n
);
1151 srenew(fr
->f_novirsum_alloc
,fr
->f_novirsum_nalloc
);
1156 fr
->f_novirsum_n
= 0;
1160 static real
cutoff_inf(real cutoff
)
1164 cutoff
= GMX_CUTOFF_INF
;
1170 static void make_adress_tf_tables(FILE *fp
,const output_env_t oenv
,
1171 t_forcerec
*fr
,const t_inputrec
*ir
,
1172 const char *tabfn
, const gmx_mtop_t
*mtop
,
1178 if (tabfn
== NULL
) {
1179 gmx_fatal(FARGS
,"No thermoforce table file given. Use -tabletf to specify a file\n");
1183 snew(fr
->atf_tabs
, ir
->adress
->n_tf_grps
);
1185 for (i
=0; i
<ir
->adress
->n_tf_grps
; i
++){
1186 j
= ir
->adress
->tf_table_index
[i
]; /* get energy group index */
1187 sprintf(buf
+ strlen(tabfn
) - strlen(ftp2ext(efXVG
)) - 1,"tf_%s.%s",
1188 *(mtop
->groups
.grpname
[mtop
->groups
.grps
[egcENER
].nm_ind
[j
]]) ,ftp2ext(efXVG
));
1189 printf("loading tf table for energygrp index %d from %s\n", ir
->adress
->tf_table_index
[j
], buf
);
1190 fr
->atf_tabs
[i
] = make_atf_table(fp
,oenv
,fr
,buf
, box
);
1195 gmx_bool
can_use_allvsall(const t_inputrec
*ir
, const gmx_mtop_t
*mtop
,
1196 gmx_bool bPrintNote
,t_commrec
*cr
,FILE *fp
)
1205 ir
->ePBC
==epbcNONE
&&
1206 ir
->vdwtype
==evdwCUT
&&
1207 ir
->coulombtype
==eelCUT
&&
1209 (ir
->implicit_solvent
== eisNO
||
1210 (ir
->implicit_solvent
==eisGBSA
&& (ir
->gb_algorithm
==egbSTILL
||
1211 ir
->gb_algorithm
==egbHCT
||
1212 ir
->gb_algorithm
==egbOBC
))) &&
1213 getenv("GMX_NO_ALLVSALL") == NULL
1216 if (bAllvsAll
&& ir
->opts
.ngener
> 1)
1218 const char *note
="NOTE: Can not use all-vs-all force loops, because there are multiple energy monitor groups; you might get significantly higher performance when using only a single energy monitor group.\n";
1224 fprintf(stderr
,"\n%s\n",note
);
1228 fprintf(fp
,"\n%s\n",note
);
1234 if(bAllvsAll
&& fp
&& MASTER(cr
))
1236 fprintf(fp
,"\nUsing accelerated all-vs-all kernels.\n\n");
1243 /* Return 1 if SSE2 support is present, otherwise 0. */
1245 forcerec_check_sse2()
1247 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE)|| defined(GMX_SSE2) )
1249 unsigned int _eax
,_ebx
,_ecx
,_edx
;
1262 #elif defined(__x86_64__)
1263 /* GCC 64-bit inline asm */
1264 __asm__ ("push %%rbx\n\tcpuid\n\tpop %%rbx\n" \
1265 : "=a" (_eax
), "=S" (_ebx
), "=c" (_ecx
), "=d" (_edx
) \
1267 #elif defined(__i386__)
1268 __asm__ ("push %%ebx\n\tcpuid\n\tpop %%ebx\n" \
1269 : "=a" (_eax
), "=S" (_ebx
), "=c" (_ecx
), "=d" (_edx
) \
1272 _eax
=_ebx
=_ecx
=_edx
=0;
1277 * SSE Bit 25 of edx should be set
1278 * SSE2 Bit 26 of edx should be set
1279 * SSE3 Bit 0 of ecx should be set
1280 * SSE4.1 Bit 19 of ecx should be set
1282 status
= (_edx
& (1 << 26)) != 0;
1287 /* Return SSE2 status */
1294 void init_forcerec(FILE *fp
,
1295 const output_env_t oenv
,
1298 const t_inputrec
*ir
,
1299 const gmx_mtop_t
*mtop
,
1300 const t_commrec
*cr
,
1307 gmx_bool bNoSolvOpt
,
1310 int i
,j
,m
,natoms
,ngrp
,negp_pp
,negptable
,egi
,egj
;
1316 gmx_bool bGenericKernelOnly
;
1317 gmx_bool bTab
,bSep14tab
,bNormalnblists
;
1319 int *nm_ind
,egp_flags
;
1321 fr
->bDomDec
= DOMAINDECOMP(cr
);
1323 natoms
= mtop
->natoms
;
1325 if (check_box(ir
->ePBC
,box
))
1327 gmx_fatal(FARGS
,check_box(ir
->ePBC
,box
));
1330 /* Test particle insertion ? */
1331 if (EI_TPI(ir
->eI
)) {
1332 /* Set to the size of the molecule to be inserted (the last one) */
1333 /* Because of old style topologies, we have to use the last cg
1334 * instead of the last molecule type.
1336 cgs
= &mtop
->moltype
[mtop
->molblock
[mtop
->nmolblock
-1].type
].cgs
;
1337 fr
->n_tpi
= cgs
->index
[cgs
->nr
] - cgs
->index
[cgs
->nr
-1];
1338 if (fr
->n_tpi
!= mtop
->mols
.index
[mtop
->mols
.nr
] - mtop
->mols
.index
[mtop
->mols
.nr
-1]) {
1339 gmx_fatal(FARGS
,"The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group.");
1345 /* Copy AdResS parameters */
1347 fr
->adress_type
= ir
->adress
->type
;
1348 fr
->adress_const_wf
= ir
->adress
->const_wf
;
1349 fr
->adress_ex_width
= ir
->adress
->ex_width
;
1350 fr
->adress_hy_width
= ir
->adress
->hy_width
;
1351 fr
->adress_icor
= ir
->adress
->icor
;
1352 fr
->adress_site
= ir
->adress
->site
;
1353 fr
->adress_ex_forcecap
= ir
->adress
->ex_forcecap
;
1354 fr
->adress_do_hybridpairs
= ir
->adress
->do_hybridpairs
;
1357 snew(fr
->adress_group_explicit
, ir
->adress
->n_energy_grps
);
1358 for (i
=0; i
< ir
->adress
->n_energy_grps
; i
++){
1359 fr
->adress_group_explicit
[i
]= ir
->adress
->group_explicit
[i
];
1362 fr
->n_adress_tf_grps
= ir
->adress
->n_tf_grps
;
1363 snew(fr
->adress_tf_table_index
, fr
->n_adress_tf_grps
);
1364 for (i
=0; i
< fr
->n_adress_tf_grps
; i
++){
1365 fr
->adress_tf_table_index
[i
]= ir
->adress
->tf_table_index
[i
];
1367 copy_rvec(ir
->adress
->refs
,fr
->adress_refs
);
1369 fr
->adress_type
= eAdressOff
;
1370 fr
->adress_do_hybridpairs
= FALSE
;
1373 /* Copy the user determined parameters */
1374 fr
->userint1
= ir
->userint1
;
1375 fr
->userint2
= ir
->userint2
;
1376 fr
->userint3
= ir
->userint3
;
1377 fr
->userint4
= ir
->userint4
;
1378 fr
->userreal1
= ir
->userreal1
;
1379 fr
->userreal2
= ir
->userreal2
;
1380 fr
->userreal3
= ir
->userreal3
;
1381 fr
->userreal4
= ir
->userreal4
;
1384 fr
->fc_stepsize
= ir
->fc_stepsize
;
1387 fr
->efep
= ir
->efep
;
1388 fr
->sc_alphavdw
= ir
->fepvals
->sc_alpha
;
1389 if (ir
->fepvals
->bScCoul
)
1391 fr
->sc_alphacoul
= ir
->fepvals
->sc_alpha
;
1392 fr
->sc_sigma6_min
= pow(ir
->fepvals
->sc_sigma_min
,6);
1396 fr
->sc_alphacoul
= 0;
1397 fr
->sc_sigma6_min
= 0; /* only needed when bScCoul is on */
1399 fr
->sc_power
= ir
->fepvals
->sc_power
;
1400 fr
->sc_r_power
= ir
->fepvals
->sc_r_power
;
1401 fr
->sc_sigma6_def
= pow(ir
->fepvals
->sc_sigma
,6);
1403 env
= getenv("GMX_SCSIGMA_MIN");
1407 sscanf(env
,"%lf",&dbl
);
1408 fr
->sc_sigma6_min
= pow(dbl
,6);
1411 fprintf(fp
,"Setting the minimum soft core sigma to %g nm\n",dbl
);
1415 bGenericKernelOnly
= FALSE
;
1416 if (getenv("GMX_NB_GENERIC") != NULL
)
1421 "Found environment variable GMX_NB_GENERIC.\n"
1422 "Disabling interaction-specific nonbonded kernels.\n\n");
1424 bGenericKernelOnly
= TRUE
;
1428 fr
->UseOptimizedKernels
= (getenv("GMX_NOOPTIMIZEDKERNELS") == NULL
);
1429 if(fp
&& fr
->UseOptimizedKernels
==FALSE
)
1432 "\nFound environment variable GMX_NOOPTIMIZEDKERNELS.\n"
1433 "Disabling SSE/SSE2/Altivec/ia64/Power6/Bluegene specific kernels.\n\n");
1436 #if ( defined(GMX_IA32_SSE2) || defined(GMX_X86_64_SSE2) || defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE)|| defined(GMX_SSE2) )
1437 if( forcerec_check_sse2() == 0 )
1439 fr
->UseOptimizedKernels
= FALSE
;
1443 /* Check if we can/should do all-vs-all kernels */
1444 fr
->bAllvsAll
= can_use_allvsall(ir
,mtop
,FALSE
,NULL
,NULL
);
1445 fr
->AllvsAll_work
= NULL
;
1446 fr
->AllvsAll_workgb
= NULL
;
1450 /* Neighbour searching stuff */
1451 fr
->bGrid
= (ir
->ns_type
== ensGRID
);
1452 fr
->ePBC
= ir
->ePBC
;
1453 fr
->bMolPBC
= ir
->bPeriodicMols
;
1454 fr
->rc_scaling
= ir
->refcoord_scaling
;
1455 copy_rvec(ir
->posres_com
,fr
->posres_com
);
1456 copy_rvec(ir
->posres_comB
,fr
->posres_comB
);
1457 fr
->rlist
= cutoff_inf(ir
->rlist
);
1458 fr
->rlistlong
= cutoff_inf(ir
->rlistlong
);
1459 fr
->eeltype
= ir
->coulombtype
;
1460 fr
->vdwtype
= ir
->vdwtype
;
1462 fr
->bTwinRange
= fr
->rlistlong
> fr
->rlist
;
1463 fr
->bEwald
= (EEL_PME(fr
->eeltype
) || fr
->eeltype
==eelEWALD
);
1465 fr
->reppow
= mtop
->ffparams
.reppow
;
1466 fr
->bvdwtab
= (fr
->vdwtype
!= evdwCUT
||
1467 !gmx_within_tol(fr
->reppow
,12.0,10*GMX_DOUBLE_EPS
));
1468 fr
->bcoultab
= (!(fr
->eeltype
== eelCUT
|| EEL_RF(fr
->eeltype
)) ||
1469 fr
->eeltype
== eelRF_ZERO
);
1471 if (getenv("GMX_FORCE_TABLES"))
1474 fr
->bcoultab
= TRUE
;
1478 fprintf(fp
,"Table routines are used for coulomb: %s\n",bool_names
[fr
->bcoultab
]);
1479 fprintf(fp
,"Table routines are used for vdw: %s\n",bool_names
[fr
->bvdwtab
]);
1482 /* Tables are used for direct ewald sum */
1485 if (EEL_PME(ir
->coulombtype
))
1488 fprintf(fp
,"Will do PME sum in reciprocal space.\n");
1489 if (ir
->coulombtype
== eelP3M_AD
)
1491 please_cite(fp
,"Hockney1988");
1492 please_cite(fp
,"Ballenegger2012");
1496 please_cite(fp
,"Essmann95a");
1499 if (ir
->ewald_geometry
== eewg3DC
)
1503 fprintf(fp
,"Using the Ewald3DC correction for systems with a slab geometry.\n");
1505 please_cite(fp
,"In-Chul99a");
1508 fr
->ewaldcoeff
=calc_ewaldcoeff(ir
->rcoulomb
, ir
->ewald_rtol
);
1509 init_ewald_tab(&(fr
->ewald_table
), cr
, ir
, fp
);
1512 fprintf(fp
,"Using a Gaussian width (1/beta) of %g nm for Ewald\n",
1517 /* Electrostatics */
1518 fr
->epsilon_r
= ir
->epsilon_r
;
1519 fr
->epsilon_rf
= ir
->epsilon_rf
;
1520 fr
->fudgeQQ
= mtop
->ffparams
.fudgeQQ
;
1521 fr
->rcoulomb_switch
= ir
->rcoulomb_switch
;
1522 fr
->rcoulomb
= cutoff_inf(ir
->rcoulomb
);
1524 /* Parameters for generalized RF */
1528 if (fr
->eeltype
== eelGRF
)
1530 init_generalized_rf(fp
,mtop
,ir
,fr
);
1532 else if (fr
->eeltype
== eelSHIFT
)
1534 for(m
=0; (m
<DIM
); m
++)
1535 box_size
[m
]=box
[m
][m
];
1537 if ((fr
->eeltype
== eelSHIFT
&& fr
->rcoulomb
> fr
->rcoulomb_switch
))
1538 set_shift_consts(fp
,fr
->rcoulomb_switch
,fr
->rcoulomb
,box_size
,fr
);
1541 fr
->bF_NoVirSum
= (EEL_FULL(fr
->eeltype
) ||
1542 gmx_mtop_ftype_count(mtop
,F_POSRES
) > 0 ||
1543 IR_ELEC_FIELD(*ir
) ||
1544 (fr
->adress_icor
!= eAdressICOff
)
1547 /* Mask that says whether or not this NBF list should be computed */
1548 /* if (fr->bMask == NULL) {
1549 ngrp = ir->opts.ngener*ir->opts.ngener;
1550 snew(fr->bMask,ngrp);*/
1551 /* Defaults to always */
1552 /* for(i=0; (i<ngrp); i++)
1553 fr->bMask[i] = TRUE;
1556 if (ncg_mtop(mtop
) > fr
->cg_nalloc
&& !DOMAINDECOMP(cr
)) {
1557 /* Count the total number of charge groups */
1558 fr
->cg_nalloc
= ncg_mtop(mtop
);
1559 srenew(fr
->cg_cm
,fr
->cg_nalloc
);
1561 if (fr
->shift_vec
== NULL
)
1562 snew(fr
->shift_vec
,SHIFTS
);
1564 if (fr
->fshift
== NULL
)
1565 snew(fr
->fshift
,SHIFTS
);
1567 if (fr
->nbfp
== NULL
) {
1568 fr
->ntype
= mtop
->ffparams
.atnr
;
1569 fr
->bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
1570 fr
->nbfp
= mk_nbfp(&mtop
->ffparams
,fr
->bBHAM
);
1573 /* Copy the energy group exclusions */
1574 fr
->egp_flags
= ir
->opts
.egp_flags
;
1576 /* Van der Waals stuff */
1577 fr
->rvdw
= cutoff_inf(ir
->rvdw
);
1578 fr
->rvdw_switch
= ir
->rvdw_switch
;
1579 if ((fr
->vdwtype
!= evdwCUT
) && (fr
->vdwtype
!= evdwUSER
) && !fr
->bBHAM
) {
1580 if (fr
->rvdw_switch
>= fr
->rvdw
)
1581 gmx_fatal(FARGS
,"rvdw_switch (%f) must be < rvdw (%f)",
1582 fr
->rvdw_switch
,fr
->rvdw
);
1584 fprintf(fp
,"Using %s Lennard-Jones, switch between %g and %g nm\n",
1585 (fr
->eeltype
==eelSWITCH
) ? "switched":"shifted",
1586 fr
->rvdw_switch
,fr
->rvdw
);
1589 if (fr
->bBHAM
&& (fr
->vdwtype
== evdwSHIFT
|| fr
->vdwtype
== evdwSWITCH
))
1590 gmx_fatal(FARGS
,"Switch/shift interaction not supported with Buckingham");
1593 fprintf(fp
,"Cut-off's: NS: %g Coulomb: %g %s: %g\n",
1594 fr
->rlist
,fr
->rcoulomb
,fr
->bBHAM
? "BHAM":"LJ",fr
->rvdw
);
1596 fr
->eDispCorr
= ir
->eDispCorr
;
1597 if (ir
->eDispCorr
!= edispcNO
)
1599 set_avcsixtwelve(fp
,fr
,mtop
);
1604 set_bham_b_max(fp
,fr
,mtop
);
1607 fr
->bGB
= (ir
->implicit_solvent
== eisGBSA
);
1608 fr
->gb_epsilon_solvent
= ir
->gb_epsilon_solvent
;
1610 /* Copy the GBSA data (radius, volume and surftens for each
1611 * atomtype) from the topology atomtype section to forcerec.
1613 snew(fr
->atype_radius
,fr
->ntype
);
1614 snew(fr
->atype_vol
,fr
->ntype
);
1615 snew(fr
->atype_surftens
,fr
->ntype
);
1616 snew(fr
->atype_gb_radius
,fr
->ntype
);
1617 snew(fr
->atype_S_hct
,fr
->ntype
);
1619 if (mtop
->atomtypes
.nr
> 0)
1621 for(i
=0;i
<fr
->ntype
;i
++)
1622 fr
->atype_radius
[i
] =mtop
->atomtypes
.radius
[i
];
1623 for(i
=0;i
<fr
->ntype
;i
++)
1624 fr
->atype_vol
[i
] = mtop
->atomtypes
.vol
[i
];
1625 for(i
=0;i
<fr
->ntype
;i
++)
1626 fr
->atype_surftens
[i
] = mtop
->atomtypes
.surftens
[i
];
1627 for(i
=0;i
<fr
->ntype
;i
++)
1628 fr
->atype_gb_radius
[i
] = mtop
->atomtypes
.gb_radius
[i
];
1629 for(i
=0;i
<fr
->ntype
;i
++)
1630 fr
->atype_S_hct
[i
] = mtop
->atomtypes
.S_hct
[i
];
1633 /* Generate the GB table if needed */
1637 fr
->gbtabscale
=2000;
1643 fr
->gbtab
=make_gb_table(fp
,oenv
,fr
,tabpfn
,fr
->gbtabscale
);
1645 init_gb(&fr
->born
,cr
,fr
,ir
,mtop
,ir
->rgbradii
,ir
->gb_algorithm
);
1647 /* Copy local gb data (for dd, this is done in dd_partition_system) */
1648 if (!DOMAINDECOMP(cr
))
1650 make_local_gb(cr
,fr
->born
,ir
->gb_algorithm
);
1654 /* Set the charge scaling */
1655 if (fr
->epsilon_r
!= 0)
1656 fr
->epsfac
= ONE_4PI_EPS0
/fr
->epsilon_r
;
1658 /* eps = 0 is infinite dieletric: no coulomb interactions */
1661 /* Reaction field constants */
1662 if (EEL_RF(fr
->eeltype
))
1663 calc_rffac(fp
,fr
->eeltype
,fr
->epsilon_r
,fr
->epsilon_rf
,
1664 fr
->rcoulomb
,fr
->temp
,fr
->zsquare
,box
,
1665 &fr
->kappa
,&fr
->k_rf
,&fr
->c_rf
);
1667 set_chargesum(fp
,fr
,mtop
);
1669 /* if we are using LR electrostatics, and they are tabulated,
1670 * the tables will contain modified coulomb interactions.
1671 * Since we want to use the non-shifted ones for 1-4
1672 * coulombic interactions, we must have an extra set of tables.
1675 /* Construct tables.
1676 * A little unnecessary to make both vdw and coul tables sometimes,
1677 * but what the heck... */
1679 bTab
= fr
->bcoultab
|| fr
->bvdwtab
;
1681 bSep14tab
= ((!bTab
|| fr
->eeltype
!=eelCUT
|| fr
->vdwtype
!=evdwCUT
||
1683 (gmx_mtop_ftype_count(mtop
,F_LJ14
) > 0 ||
1684 gmx_mtop_ftype_count(mtop
,F_LJC14_Q
) > 0 ||
1685 gmx_mtop_ftype_count(mtop
,F_LJC_PAIRS_NB
) > 0));
1687 negp_pp
= ir
->opts
.ngener
- ir
->nwall
;
1690 bNormalnblists
= TRUE
;
1693 bNormalnblists
= (ir
->eDispCorr
!= edispcNO
);
1694 for(egi
=0; egi
<negp_pp
; egi
++) {
1695 for(egj
=egi
; egj
<negp_pp
; egj
++) {
1696 egp_flags
= ir
->opts
.egp_flags
[GID(egi
,egj
,ir
->opts
.ngener
)];
1697 if (!(egp_flags
& EGP_EXCL
)) {
1698 if (egp_flags
& EGP_TABLE
) {
1701 bNormalnblists
= TRUE
;
1706 if (bNormalnblists
) {
1707 fr
->nnblists
= negptable
+ 1;
1709 fr
->nnblists
= negptable
;
1711 if (fr
->nnblists
> 1)
1712 snew(fr
->gid2nblists
,ir
->opts
.ngener
*ir
->opts
.ngener
);
1714 snew(fr
->nblists
,fr
->nnblists
);
1716 /* This code automatically gives table length tabext without cut-off's,
1717 * in that case grompp should already have checked that we do not need
1718 * normal tables and we only generate tables for 1-4 interactions.
1720 rtab
= ir
->rlistlong
+ ir
->tabext
;
1723 /* make tables for ordinary interactions */
1724 if (bNormalnblists
) {
1725 make_nbf_tables(fp
,oenv
,fr
,rtab
,cr
,tabfn
,NULL
,NULL
,&fr
->nblists
[0]);
1727 fr
->tab14
= fr
->nblists
[0].tab
;
1732 if (negptable
> 0) {
1733 /* Read the special tables for certain energy group pairs */
1734 nm_ind
= mtop
->groups
.grps
[egcENER
].nm_ind
;
1735 for(egi
=0; egi
<negp_pp
; egi
++) {
1736 for(egj
=egi
; egj
<negp_pp
; egj
++) {
1737 egp_flags
= ir
->opts
.egp_flags
[GID(egi
,egj
,ir
->opts
.ngener
)];
1738 if ((egp_flags
& EGP_TABLE
) && !(egp_flags
& EGP_EXCL
)) {
1739 nbl
= &(fr
->nblists
[m
]);
1740 if (fr
->nnblists
> 1) {
1741 fr
->gid2nblists
[GID(egi
,egj
,ir
->opts
.ngener
)] = m
;
1743 /* Read the table file with the two energy groups names appended */
1744 make_nbf_tables(fp
,oenv
,fr
,rtab
,cr
,tabfn
,
1745 *mtop
->groups
.grpname
[nm_ind
[egi
]],
1746 *mtop
->groups
.grpname
[nm_ind
[egj
]],
1749 } else if (fr
->nnblists
> 1) {
1750 fr
->gid2nblists
[GID(egi
,egj
,ir
->opts
.ngener
)] = 0;
1758 /* generate extra tables with plain Coulomb for 1-4 interactions only */
1759 fr
->tab14
= make_tables(fp
,oenv
,fr
,MASTER(cr
),tabpfn
,rtab
,
1760 GMX_MAKETABLES_14ONLY
);
1763 /* Read AdResS Thermo Force table if needed */
1764 if(fr
->adress_icor
== eAdressICThermoForce
)
1766 /* old todo replace */
1768 if (ir
->adress
->n_tf_grps
> 0){
1769 make_adress_tf_tables(fp
,oenv
,fr
,ir
,tabfn
, mtop
, box
);
1772 /* load the default table */
1773 snew(fr
->atf_tabs
, 1);
1774 fr
->atf_tabs
[DEFAULT_TF_TABLE
] = make_atf_table(fp
,oenv
,fr
,tabafn
, box
);
1779 fr
->nwall
= ir
->nwall
;
1780 if (ir
->nwall
&& ir
->wall_type
==ewtTABLE
)
1782 make_wall_tables(fp
,oenv
,ir
,tabfn
,&mtop
->groups
,fr
);
1785 if (fcd
&& tabbfn
) {
1786 fcd
->bondtab
= make_bonded_tables(fp
,
1787 F_TABBONDS
,F_TABBONDSNC
,
1789 fcd
->angletab
= make_bonded_tables(fp
,
1792 fcd
->dihtab
= make_bonded_tables(fp
,
1797 fprintf(debug
,"No fcdata or table file name passed, can not read table, can not do bonded interactions\n");
1800 if (ir
->eDispCorr
!= edispcNO
)
1802 calc_enervirdiff(fp
,ir
->eDispCorr
,fr
);
1805 /* QM/MM initialization if requested
1809 fprintf(stderr
,"QM/MM calculation requested.\n");
1812 fr
->bQMMM
= ir
->bQMMM
;
1813 fr
->qr
= mk_QMMMrec();
1815 /* Set all the static charge group info */
1816 fr
->cginfo_mb
= init_cginfo_mb(fp
,mtop
,fr
,bNoSolvOpt
);
1818 if (DOMAINDECOMP(cr
)) {
1821 fr
->cginfo
= cginfo_expand(mtop
->nmolblock
,fr
->cginfo_mb
);
1824 if (!DOMAINDECOMP(cr
))
1826 /* When using particle decomposition, the effect of the second argument,
1827 * which sets fr->hcg, is corrected later in do_md and init_em.
1829 forcerec_set_ranges(fr
,ncg_mtop(mtop
),ncg_mtop(mtop
),
1830 mtop
->natoms
,mtop
->natoms
,mtop
->natoms
);
1833 fr
->print_force
= print_force
;
1836 /* coarse load balancing vars */
1841 /* Initialize neighbor search */
1842 init_ns(fp
,cr
,&fr
->ns
,fr
,mtop
,box
);
1844 if (cr
->duty
& DUTY_PP
){
1845 gmx_setup_kernels(fp
,bGenericKernelOnly
);
1847 gmx_setup_adress_kernels(fp
,bGenericKernelOnly
);
1851 #define pr_real(fp,r) fprintf(fp,"%s: %e\n",#r,r)
1852 #define pr_int(fp,i) fprintf((fp),"%s: %d\n",#i,i)
1853 #define pr_bool(fp,b) fprintf((fp),"%s: %s\n",#b,bool_names[b])
1855 void pr_forcerec(FILE *fp
,t_forcerec
*fr
,t_commrec
*cr
)
1859 pr_real(fp
,fr
->rlist
);
1860 pr_real(fp
,fr
->rcoulomb
);
1861 pr_real(fp
,fr
->fudgeQQ
);
1862 pr_bool(fp
,fr
->bGrid
);
1863 pr_bool(fp
,fr
->bTwinRange
);
1864 /*pr_int(fp,fr->cg0);
1865 pr_int(fp,fr->hcg);*/
1866 for(i
=0; i
<fr
->nnblists
; i
++)
1867 pr_int(fp
,fr
->nblists
[i
].tab
.n
);
1868 pr_real(fp
,fr
->rcoulomb_switch
);
1869 pr_real(fp
,fr
->rcoulomb
);