2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
5 * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "gromacs/math/utilities.h"
41 #include "gromacs/pbcutil/pbc.h"
42 #include "gromacs/legacyheaders/types/simple.h"
43 #include "gromacs/legacyheaders/typedefs.h"
44 #include "gromacs/math/vec.h"
46 #include "gromacs/utility/fatalerror.h"
58 real l2 = adressr+adressw;
67 pbc_dx(pbc, (*ref), x, dx);
71 rvec_sub((*ref), x, dx);
77 /* default to explicit simulation */
80 /* constant value for weighting function = adressw */
81 return fr->adress_const_wf;
83 /* plane through center of ref, varies in x direction */
87 /* point at center of ref, assuming cubic geometry */
88 for (i = 0; i < 3; i++)
90 sqr_dl += dx[i]*dx[i];
94 /* default to explicit simulation */
100 /* molecule is coarse grained */
105 /* molecule is explicit */
106 else if (dl < adressr)
113 tmp = cos((dl-adressr)*M_PI/2/adressw);
119 update_adress_weights_com(FILE gmx_unused * fplog,
128 int icg, k, k0, k1, d;
129 real nrcg, inv_ncg, mtot, inv_mtot;
133 real adressr, adressw;
139 int n_hyb, n_ex, n_cg;
145 adresstype = fr->adress_type;
146 adressr = fr->adress_ex_width;
147 adressw = fr->adress_hy_width;
148 massT = mdatoms->massT;
150 ref = &(fr->adress_refs);
153 /* Since this is center of mass AdResS, the vsite is not guaranteed
154 * to be on the same node as the constructing atoms. Therefore we
155 * loop over the charge groups, calculate their center of mass,
156 * then use this to calculate wf for each atom. This wastes vsite
157 * construction, but it's the only way to assure that the explicit
158 * atoms have the same wf as their vsite. */
161 fprintf(fplog, "Calculating center of mass for charge groups %d to %d\n",
164 cgindex = cgs->index;
166 /* Compute the center of mass for all charge groups */
167 for (icg = cg0; (icg < cg1); icg++)
174 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
179 else if (wf[k0] == 1)
191 for (k = k0; (k < k1); k++)
200 for (k = k0; (k < k1); k++)
202 for (d = 0; (d < DIM); d++)
204 ix[d] += x[k][d]*massT[k];
207 for (d = 0; (d < DIM); d++)
212 /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
218 for (k = k0; (k < k1); k++)
220 for (d = 0; (d < DIM); d++)
225 for (d = 0; (d < DIM); d++)
231 /* Set wf of all atoms in charge group equal to wf of com */
232 wf[k0] = adress_weight(ix, adresstype, adressr, adressw, ref, pbc, fr);
238 else if (wf[k0] == 1)
247 for (k = (k0+1); (k < k1); k++)
255 void update_adress_weights_atom_per_atom(
264 int icg, k, k0, k1, d;
265 real nrcg, inv_ncg, mtot, inv_mtot;
269 real adressr, adressw;
275 int n_hyb, n_ex, n_cg;
281 adresstype = fr->adress_type;
282 adressr = fr->adress_ex_width;
283 adressw = fr->adress_hy_width;
284 massT = mdatoms->massT;
286 ref = &(fr->adress_refs);
288 cgindex = cgs->index;
290 /* Weighting function is determined for each atom individually.
291 * This is an approximation
292 * as in the theory requires an interpolation based on the center of masses.
293 * Should be used with caution */
295 for (icg = cg0; (icg < cg1); icg++)
298 k1 = cgindex[icg + 1];
301 for (k = (k0); (k < k1); k++)
303 wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
322 update_adress_weights_cog(t_iparams ip[],
329 int i, j, k, nr, nra, inc;
330 int ftype, adresstype;
331 t_iatom avsite, ai, aj, ak, al;
333 real adressr, adressw;
336 int n_hyb, n_ex, n_cg;
338 adresstype = fr->adress_type;
339 adressr = fr->adress_ex_width;
340 adressw = fr->adress_hy_width;
342 ref = &(fr->adress_refs);
350 /* Since this is center of geometry AdResS, we know the vsite
351 * is in the same charge group node as the constructing atoms.
352 * Loop over vsite types, calculate the weight of the vsite,
353 * then assign that weight to the constructing atoms. */
355 for (ftype = 0; (ftype < F_NRE); ftype++)
357 if (interaction_function[ftype].flags & IF_VSITE)
359 nra = interaction_function[ftype].nratoms;
360 nr = ilist[ftype].nr;
361 ia = ilist[ftype].iatoms;
363 for (i = 0; (i < nr); )
365 /* The vsite and first constructing atom */
368 wf[avsite] = adress_weight(x[avsite], adresstype, adressr, adressw, ref, pbc, fr);
375 else if (wf[ai] == 1)
384 /* Assign the vsite wf to rest of constructing atoms depending on type */
433 inc = 3*ip[ia[0]].vsiten.n;
434 for (j = 3; j < inc; j += 3)
441 gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
442 ftype, __FILE__, __LINE__);
445 /* Increment loop variables */
454 update_adress_weights_atom(int cg0,
465 real adressr, adressw;
470 adresstype = fr->adress_type;
471 adressr = fr->adress_ex_width;
472 adressw = fr->adress_hy_width;
473 massT = mdatoms->massT;
475 ref = &(fr->adress_refs);
476 cgindex = cgs->index;
478 /* Only use first atom in charge group.
479 * We still can't be sure that the vsite and constructing
480 * atoms are on the same processor, so we must calculate
481 * in the same way as com adress. */
483 for (icg = cg0; (icg < cg1); icg++)
487 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
489 /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
490 for (k = (k0+1); (k < k1); k++)
498 adress_thermo_force(int start,
507 int iatom, n0, nnn, nrcg, i;
509 real adressw, adressr;
511 unsigned short * ptype;
517 real w, wsq, wmin1, wmin1sq, wp, wt, rinv, sqr_dl, dl;
518 real eps, eps2, F, Geps, Heps2, Fp, dmu_dwp, dwp_dr, fscal;
520 adresstype = fr->adress_type;
521 adressw = fr->adress_hy_width;
522 adressr = fr->adress_ex_width;
523 cgindex = cgs->index;
524 ptype = mdatoms->ptype;
525 ref = &(fr->adress_refs);
528 for (iatom = start; (iatom < start+homenr); iatom++)
530 if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
532 if (ptype[iatom] == eptVSite)
535 /* is it hybrid or apply the thermodynamics force everywhere?*/
536 if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
538 if (fr->n_adress_tf_grps > 0)
540 /* multi component tf is on, select the right table */
541 ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
542 tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
546 /* just on component*/
547 ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].data;
548 tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
554 pbc_dx(pbc, (*ref), x[iatom], dr);
558 rvec_sub((*ref), x[iatom], dr);
564 /* calculate distace to adress center again */
569 /* plane through center of ref, varies in x direction */
570 sqr_dl = dr[0]*dr[0];
571 rinv = gmx_invsqrt(dr[0]*dr[0]);
574 /* point at center of ref, assuming cubic geometry */
575 for (i = 0; i < 3; i++)
577 sqr_dl += dr[i]*dr[i];
579 rinv = gmx_invsqrt(sqr_dl);
582 /* This case should not happen */
587 /* table origin is adress center */
594 Geps = eps*ATFtab[nnn+2];
595 Heps2 = eps2*ATFtab[nnn+3];
597 F = (Fp+Geps+2.0*Heps2)*tabscale;
601 f[iatom][0] += fscal*dr[0];
602 if (adresstype != eAdressXSplit)
604 f[iatom][1] += fscal*dr[1];
605 f[iatom][2] += fscal*dr[2];
613 gmx_bool egp_explicit(t_forcerec * fr, int egp_nr)
615 return fr->adress_group_explicit[egp_nr];
618 gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr)
620 return !fr->adress_group_explicit[egp_nr];