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 */
83 for (i = 0; i < 3; i++)
85 sqr_dl += dx[i]*dx[i];
89 /* default to explicit simulation */
95 /* molecule is coarse grained */
100 /* molecule is explicit */
101 else if (dl < adressr)
108 tmp = cos((dl-adressr)*M_PI/2/adressw);
114 update_adress_weights_com(FILE gmx_unused * fplog,
123 int icg, k, k0, k1, d;
124 real nrcg, inv_ncg, mtot, inv_mtot;
128 real adressr, adressw;
134 int n_hyb, n_ex, n_cg;
140 adresstype = fr->adress_type;
141 adressr = fr->adress_ex_width;
142 adressw = fr->adress_hy_width;
143 massT = mdatoms->massT;
145 ref = &(fr->adress_refs);
148 /* Since this is center of mass AdResS, the vsite is not guaranteed
149 * to be on the same node as the constructing atoms. Therefore we
150 * loop over the charge groups, calculate their center of mass,
151 * then use this to calculate wf for each atom. This wastes vsite
152 * construction, but it's the only way to assure that the explicit
153 * atoms have the same wf as their vsite. */
156 fprintf(fplog, "Calculating center of mass for charge groups %d to %d\n",
159 cgindex = cgs->index;
161 /* Compute the center of mass for all charge groups */
162 for (icg = cg0; (icg < cg1); icg++)
169 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
174 else if (wf[k0] == 1)
186 for (k = k0; (k < k1); k++)
195 for (k = k0; (k < k1); k++)
197 for (d = 0; (d < DIM); d++)
199 ix[d] += x[k][d]*massT[k];
202 for (d = 0; (d < DIM); d++)
207 /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
213 for (k = k0; (k < k1); k++)
215 for (d = 0; (d < DIM); d++)
220 for (d = 0; (d < DIM); d++)
226 /* Set wf of all atoms in charge group equal to wf of com */
227 wf[k0] = adress_weight(ix, adresstype, adressr, adressw, ref, pbc, fr);
233 else if (wf[k0] == 1)
242 for (k = (k0+1); (k < k1); k++)
250 void update_adress_weights_atom_per_atom(
259 int icg, k, k0, k1, d;
260 real nrcg, inv_ncg, mtot, inv_mtot;
264 real adressr, adressw;
270 int n_hyb, n_ex, n_cg;
276 adresstype = fr->adress_type;
277 adressr = fr->adress_ex_width;
278 adressw = fr->adress_hy_width;
279 massT = mdatoms->massT;
281 ref = &(fr->adress_refs);
283 cgindex = cgs->index;
285 /* Weighting function is determined for each atom individually.
286 * This is an approximation
287 * as in the theory requires an interpolation based on the center of masses.
288 * Should be used with caution */
290 for (icg = cg0; (icg < cg1); icg++)
293 k1 = cgindex[icg + 1];
296 for (k = (k0); (k < k1); k++)
298 wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
317 update_adress_weights_cog(t_iparams ip[],
324 int i, j, k, nr, nra, inc;
325 int ftype, adresstype;
326 t_iatom avsite, ai, aj, ak, al;
328 real adressr, adressw;
331 int n_hyb, n_ex, n_cg;
333 adresstype = fr->adress_type;
334 adressr = fr->adress_ex_width;
335 adressw = fr->adress_hy_width;
337 ref = &(fr->adress_refs);
345 /* Since this is center of geometry AdResS, we know the vsite
346 * is in the same charge group node as the constructing atoms.
347 * Loop over vsite types, calculate the weight of the vsite,
348 * then assign that weight to the constructing atoms. */
350 for (ftype = 0; (ftype < F_NRE); ftype++)
352 if (interaction_function[ftype].flags & IF_VSITE)
354 nra = interaction_function[ftype].nratoms;
355 nr = ilist[ftype].nr;
356 ia = ilist[ftype].iatoms;
358 for (i = 0; (i < nr); )
360 /* The vsite and first constructing atom */
363 wf[avsite] = adress_weight(x[avsite], adresstype, adressr, adressw, ref, pbc, fr);
370 else if (wf[ai] == 1)
379 /* Assign the vsite wf to rest of constructing atoms depending on type */
428 inc = 3*ip[ia[0]].vsiten.n;
429 for (j = 3; j < inc; j += 3)
436 gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
437 ftype, __FILE__, __LINE__);
440 /* Increment loop variables */
449 update_adress_weights_atom(int cg0,
460 real adressr, adressw;
465 adresstype = fr->adress_type;
466 adressr = fr->adress_ex_width;
467 adressw = fr->adress_hy_width;
468 massT = mdatoms->massT;
470 ref = &(fr->adress_refs);
471 cgindex = cgs->index;
473 /* Only use first atom in charge group.
474 * We still can't be sure that the vsite and constructing
475 * atoms are on the same processor, so we must calculate
476 * in the same way as com adress. */
478 for (icg = cg0; (icg < cg1); icg++)
482 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
484 /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
485 for (k = (k0+1); (k < k1); k++)
493 adress_thermo_force(int start,
502 int iatom, n0, nnn, nrcg, i;
504 real adressw, adressr;
506 unsigned short * ptype;
512 real w, wsq, wmin1, wmin1sq, wp, wt, rinv, sqr_dl, dl;
513 real eps, eps2, F, Geps, Heps2, Fp, dmu_dwp, dwp_dr, fscal;
515 adresstype = fr->adress_type;
516 adressw = fr->adress_hy_width;
517 adressr = fr->adress_ex_width;
518 cgindex = cgs->index;
519 ptype = mdatoms->ptype;
520 ref = &(fr->adress_refs);
523 for (iatom = start; (iatom < start+homenr); iatom++)
525 if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
527 if (ptype[iatom] == eptVSite)
530 /* is it hybrid or apply the thermodynamics force everywhere?*/
531 if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
533 if (fr->n_adress_tf_grps > 0)
535 /* multi component tf is on, select the right table */
536 ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
537 tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
541 /* just on component*/
542 ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].data;
543 tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
549 pbc_dx(pbc, (*ref), x[iatom], dr);
553 rvec_sub((*ref), x[iatom], dr);
559 /* calculate distace to adress center again */
564 /* plane through center of ref, varies in x direction */
565 sqr_dl = dr[0]*dr[0];
566 rinv = gmx_invsqrt(dr[0]*dr[0]);
569 /* point at center of ref, assuming cubic geometry */
570 for (i = 0; i < 3; i++)
572 sqr_dl += dr[i]*dr[i];
574 rinv = gmx_invsqrt(sqr_dl);
577 /* This case should not happen */
582 /* table origin is adress center */
589 Geps = eps*ATFtab[nnn+2];
590 Heps2 = eps2*ATFtab[nnn+3];
592 F = (Fp+Geps+2.0*Heps2)*tabscale;
596 f[iatom][0] += fscal*dr[0];
597 if (adresstype != eAdressXSplit)
599 f[iatom][1] += fscal*dr[1];
600 f[iatom][2] += fscal*dr[2];
608 gmx_bool egp_explicit(t_forcerec * fr, int egp_nr)
610 return fr->adress_group_explicit[egp_nr];
613 gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr)
615 return !fr->adress_group_explicit[egp_nr];