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 Christoph Junghans, Brad Lambeth, and possibly others.
12 * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
13 * All rights reserved.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
39 #include "types/simple.h"
53 real l2 = adressr+adressw;
62 pbc_dx(pbc,(*ref),x,dx);
66 rvec_sub((*ref),x,dx);
72 /* default to explicit simulation */
75 /* constant value for weighting function = adressw */
76 return fr->adress_const_wf;
78 /* plane through center of ref, varies in x direction */
82 /* point at center of ref, assuming cubic geometry */
84 sqr_dl += dx[i]*dx[i];
88 /* default to explicit simulation */
94 /* molecule is coarse grained */
99 /* molecule is explicit */
100 else if (dl < adressr)
107 tmp=cos((dl-adressr)*M_PI/2/adressw);
113 update_adress_weights_com(FILE * fplog,
123 real nrcg,inv_ncg,mtot,inv_mtot;
127 real adressr,adressw;
133 int n_hyb, n_ex, n_cg;
139 adresstype = fr->adress_type;
140 adressr = fr->adress_ex_width;
141 adressw = fr->adress_hy_width;
142 massT = mdatoms->massT;
144 ref = &(fr->adress_refs);
147 /* Since this is center of mass AdResS, the vsite is not guaranteed
148 * to be on the same node as the constructing atoms. Therefore we
149 * loop over the charge groups, calculate their center of mass,
150 * then use this to calculate wf for each atom. This wastes vsite
151 * construction, but it's the only way to assure that the explicit
152 * atoms have the same wf as their vsite. */
155 fprintf(fplog,"Calculating center of mass for charge groups %d to %d\n",
158 cgindex = cgs->index;
160 /* Compute the center of mass for all charge groups */
161 for(icg=cg0; (icg<cg1); icg++)
168 wf[k0] = adress_weight(x[k0],adresstype,adressr,adressw,ref,pbc,fr);
169 if (wf[k0]==0){ n_cg++;}
170 else if (wf[k0]==1){ n_ex++;}
176 for(k=k0; (k<k1); k++)
185 for(k=k0; (k<k1); k++)
187 for(d=0; (d<DIM); d++)
189 ix[d] += x[k][d]*massT[k];
192 for(d=0; (d<DIM); d++)
197 /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
203 for(k=k0; (k<k1); k++)
205 for(d=0; (d<DIM); d++)
210 for(d=0; (d<DIM); d++)
216 /* Set wf of all atoms in charge group equal to wf of com */
217 wf[k0] = adress_weight(ix,adresstype,adressr,adressw,ref,pbc, fr);
219 if (wf[k0]==0){ n_cg++;}
220 else if (wf[k0]==1){ n_ex++;}
223 for(k=(k0+1); (k<k1); k++)
231 adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
235 void update_adress_weights_atom_per_atom(
245 real nrcg,inv_ncg,mtot,inv_mtot;
249 real adressr,adressw;
255 int n_hyb, n_ex, n_cg;
261 adresstype = fr->adress_type;
262 adressr = fr->adress_ex_width;
263 adressw = fr->adress_hy_width;
264 massT = mdatoms->massT;
266 ref = &(fr->adress_refs);
268 cgindex = cgs->index;
270 /* Weighting function is determined for each atom individually.
271 * This is an approximation
272 * as in the theory requires an interpolation based on the center of masses.
273 * Should be used with caution */
275 for (icg = cg0; (icg < cg1); icg++) {
277 k1 = cgindex[icg + 1];
280 for (k = (k0); (k < k1); k++) {
281 wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
284 } else if (wf[k] == 1) {
292 adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
296 update_adress_weights_cog(t_iparams ip[],
303 int i,j,k,nr,nra,inc;
304 int ftype,adresstype;
305 t_iatom avsite,ai,aj,ak,al;
307 real adressr,adressw;
310 int n_hyb, n_ex, n_cg;
312 adresstype = fr->adress_type;
313 adressr = fr->adress_ex_width;
314 adressw = fr->adress_hy_width;
316 ref = &(fr->adress_refs);
324 /* Since this is center of geometry AdResS, we know the vsite
325 * is in the same charge group node as the constructing atoms.
326 * Loop over vsite types, calculate the weight of the vsite,
327 * then assign that weight to the constructing atoms. */
329 for(ftype=0; (ftype<F_NRE); ftype++)
331 if (interaction_function[ftype].flags & IF_VSITE)
333 nra = interaction_function[ftype].nratoms;
334 nr = ilist[ftype].nr;
335 ia = ilist[ftype].iatoms;
339 /* The vsite and first constructing atom */
342 wf[avsite] = adress_weight(x[avsite],adresstype,adressr,adressw,ref,pbc,fr);
347 } else if (wf[ai] == 1) {
353 /* Assign the vsite wf to rest of constructing atoms depending on type */
401 inc = 3*ip[ia[0]].vsiten.n;
402 for(j=3; j<inc; j+=3)
409 gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
410 ftype,__FILE__,__LINE__);
413 /* Increment loop variables */
420 adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
424 update_adress_weights_atom(int cg0,
435 real adressr,adressw;
440 adresstype = fr->adress_type;
441 adressr = fr->adress_ex_width;
442 adressw = fr->adress_hy_width;
443 massT = mdatoms->massT;
445 ref = &(fr->adress_refs);
446 cgindex = cgs->index;
448 /* Only use first atom in charge group.
449 * We still can't be sure that the vsite and constructing
450 * atoms are on the same processor, so we must calculate
451 * in the same way as com adress. */
453 for(icg=cg0; (icg<cg1); icg++)
457 wf[k0] = adress_weight(x[k0],adresstype,adressr,adressw,ref,pbc,fr);
459 /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
460 for(k=(k0+1); (k<k1); k++)
467 void adress_set_kernel_flags(int n_ex, int n_hyb, int n_cg, t_mdatoms * mdatoms){
469 /* With domain decomposition we can check weather a cpu calculates only
470 * coarse-grained or explicit interactions. If so we use standard gromacs kernels
471 * on this proc. See also nonbonded.c */
473 if (n_hyb ==0 && n_ex == 0){
474 /* all particles on this proc are coarse-grained, use standard gromacs kernels */
475 if (!mdatoms->purecg){
476 mdatoms->purecg = TRUE;
477 if (debug) fprintf (debug, "adress.c: pure cg kernels on this proc\n");
482 if (mdatoms->purecg){
483 /* now this processor has hybrid particles again, call the hybrid kernels */
484 mdatoms->purecg = FALSE;
488 if (n_hyb ==0 && n_cg == 0){
489 /* all particles on this proc are atomistic, use standard gromacs kernels */
490 if (!mdatoms->pureex){
491 mdatoms->pureex = TRUE;
492 if (debug) fprintf (debug, "adress.c: pure ex kernels on this proc\n");
497 if (mdatoms->pureex){
498 mdatoms->pureex = FALSE;
504 adress_thermo_force(int start,
513 int iatom,n0,nnn,nrcg, i;
515 real adressw, adressr;
517 unsigned short * ptype;
523 real w,wsq,wmin1,wmin1sq,wp,wt,rinv, sqr_dl, dl;
524 real eps,eps2,F,Geps,Heps2,Fp,dmu_dwp,dwp_dr,fscal;
526 adresstype = fr->adress_type;
527 adressw = fr->adress_hy_width;
528 adressr = fr->adress_ex_width;
529 cgindex = cgs->index;
530 ptype = mdatoms->ptype;
531 ref = &(fr->adress_refs);
534 for(iatom=start; (iatom<start+homenr); iatom++)
536 if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
538 if (ptype[iatom] == eptVSite)
541 /* is it hybrid or apply the thermodynamics force everywhere?*/
542 if ( mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
544 if (fr->n_adress_tf_grps > 0 ){
545 /* multi component tf is on, select the right table */
546 ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].tab;
547 tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
550 /* just on component*/
551 ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].tab;
552 tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
558 pbc_dx(pbc,(*ref),x[iatom],dr);
562 rvec_sub((*ref),x[iatom],dr);
568 /* calculate distace to adress center again */
573 /* plane through center of ref, varies in x direction */
574 sqr_dl = dr[0]*dr[0];
575 rinv = gmx_invsqrt(dr[0]*dr[0]);
578 /* point at center of ref, assuming cubic geometry */
580 sqr_dl += dr[i]*dr[i];
582 rinv = gmx_invsqrt(sqr_dl);
585 /* This case should not happen */
590 /* table origin is adress center */
597 Geps = eps*ATFtab[nnn+2];
598 Heps2 = eps2*ATFtab[nnn+3];
600 F = (Fp+Geps+2.0*Heps2)*tabscale;
604 f[iatom][0] += fscal*dr[0];
605 if (adresstype != eAdressXSplit)
607 f[iatom][1] += fscal*dr[1];
608 f[iatom][2] += fscal*dr[2];
616 gmx_bool egp_explicit(t_forcerec * fr, int egp_nr)
618 return fr->adress_group_explicit[egp_nr];
621 gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr)
623 return !fr->adress_group_explicit[egp_nr];