2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/fileio/filetypes.h"
45 #include "gromacs/gmxlib/nrnb.h"
46 #include "gromacs/math/utilities.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/mdlib/force.h"
49 #include "gromacs/mdtypes/inputrec.h"
50 #include "gromacs/mdtypes/md_enums.h"
51 #include "gromacs/mdtypes/mdatom.h"
52 #include "gromacs/mdtypes/nblist.h"
53 #include "gromacs/tables/forcetable.h"
54 #include "gromacs/topology/topology.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
59 void make_wall_tables(FILE *fplog
,
60 const t_inputrec
*ir
, const char *tabfn
,
61 const gmx_groups_t
*groups
,
68 negp_pp
= ir
->opts
.ngener
- ir
->nwall
;
69 nm_ind
= groups
->grps
[egcENER
].nm_ind
;
73 fprintf(fplog
, "Reading user tables for %d energy groups with %d walls\n",
77 snew(fr
->wall_tab
, ir
->nwall
);
78 for (int w
= 0; w
< ir
->nwall
; w
++)
80 snew(fr
->wall_tab
[w
], negp_pp
);
81 for (int egp
= 0; egp
< negp_pp
; egp
++)
83 /* If the energy group pair is excluded, we don't need a table */
84 if (!(fr
->egp_flags
[egp
*ir
->opts
.ngener
+negp_pp
+w
] & EGP_EXCL
))
86 sprintf(buf
, "%s", tabfn
);
87 sprintf(buf
+ strlen(tabfn
) - strlen(ftp2ext(efXVG
)) - 1, "_%s_%s.%s",
88 *groups
->grpname
[nm_ind
[egp
]],
89 *groups
->grpname
[nm_ind
[negp_pp
+w
]],
91 fr
->wall_tab
[w
][egp
] = make_tables(fplog
, fr
->ic
, buf
, 0,
92 GMX_MAKETABLES_FORCEUSER
);
94 /* Since wall have no charge, we can compress the table */
95 for (int i
= 0; i
<= fr
->wall_tab
[w
][egp
]->n
; i
++)
97 for (int j
= 0; j
< 8; j
++)
99 fr
->wall_tab
[w
][egp
]->data
[8*i
+j
] =
100 fr
->wall_tab
[w
][egp
]->data
[12*i
+4+j
];
108 static void wall_error(int a
, const rvec
*x
, real r
)
111 "An atom is beyond the wall: coordinates %f %f %f, distance %f\n"
112 "You might want to use the mdp option wall_r_linpot",
113 x
[a
][XX
], x
[a
][YY
], x
[a
][ZZ
], r
);
116 real
do_walls(const t_inputrec
*ir
, t_forcerec
*fr
, matrix box
, const t_mdatoms
*md
,
117 const rvec x
[], rvec f
[], real lambda
, real Vlj
[], t_nrnb
*nrnb
)
120 int ntw
[2], at
, ntype
, ngid
, ggid
, *egp_flags
, *type
;
121 real
*nbfp
, lamfac
, fac_d
[2], fac_r
[2], Cd
, Cr
;
122 real wall_z
[2], r
, mr
, r1
, r2
, r4
, Vd
, Vr
, V
= 0, Fd
, Fr
, F
= 0, dvdlambda
;
125 real tabscale
, *VFtab
, rt
, eps
, eps2
, Yt
, Ft
, Geps
, Heps2
, Fp
, VV
, FF
;
126 unsigned short *gid
= md
->cENER
;
130 ngid
= ir
->opts
.ngener
;
133 egp_flags
= fr
->egp_flags
;
135 for (int w
= 0; w
< nwall
; w
++)
137 ntw
[w
] = 2*ntype
*ir
->wall_atomtype
[w
];
138 switch (ir
->wall_type
)
141 fac_d
[w
] = ir
->wall_density
[w
]*M_PI
/6;
142 fac_r
[w
] = ir
->wall_density
[w
]*M_PI
/45;
145 fac_d
[w
] = ir
->wall_density
[w
]*M_PI
/2;
146 fac_r
[w
] = ir
->wall_density
[w
]*M_PI
/5;
153 wall_z
[1] = box
[ZZ
][ZZ
];
157 for (int lam
= 0; lam
< (md
->nPerturbed
? 2 : 1); lam
++)
179 for (int i
= 0; i
< md
->homenr
; i
++)
181 for (int w
= 0; w
< std::min(nwall
, 2); w
++)
183 /* The wall energy groups are always at the end of the list */
184 ggid
= gid
[i
]*ngid
+ ngid
- nwall
+ w
;
186 /* nbfp now includes the 6/12 derivative prefactors */
187 Cd
= nbfp
[ntw
[w
]+2*at
]/6;
188 Cr
= nbfp
[ntw
[w
]+2*at
+1]/12;
189 if (!((Cd
== 0 && Cr
== 0) || (egp_flags
[ggid
] & EGP_EXCL
)))
197 r
= wall_z
[1] - x
[i
][ZZ
];
199 if (r
< ir
->wall_r_linpot
)
201 mr
= ir
->wall_r_linpot
- r
;
202 r
= ir
->wall_r_linpot
;
208 switch (ir
->wall_type
)
215 tab
= fr
->wall_tab
[w
][gid
[i
]];
216 tabscale
= tab
->scale
;
220 n0
= static_cast<int>(rt
);
223 /* Beyond the table range, set V and F to zero */
235 Geps
= VFtab
[nnn
+2]*eps
;
236 Heps2
= VFtab
[nnn
+3]*eps2
;
237 Fp
= Ft
+ Geps
+ Heps2
;
239 FF
= Fp
+ Geps
+ 2.0*Heps2
;
246 Geps
= VFtab
[nnn
+2]*eps
;
247 Heps2
= VFtab
[nnn
+3]*eps2
;
248 Fp
= Ft
+ Geps
+ Heps2
;
250 FF
= Fp
+ Geps
+ 2.0*Heps2
;
254 F
= -lamfac
*(Fd
+ Fr
)*tabscale
;
265 Vd
= fac_d
[w
]*Cd
*r2
*r1
;
266 Vr
= fac_r
[w
]*Cr
*r4
*r4
*r1
;
268 F
= lamfac
*(9*Vr
- 3*Vd
)*r1
;
279 Vr
= fac_r
[w
]*Cr
*r4
*r4
*r2
;
281 F
= lamfac
*(10*Vr
- 4*Vd
)*r1
;
294 F
= lamfac
*(12*Vr
- 6*Vd
)*r1
;
307 Vlj
[ggid
] += lamfac
*V
;
310 /* Because of the single sum virial calculation we need
311 * to add the full virial contribution of the walls.
312 * Since the force only has a z-component, there is only
313 * a contribution to the z component of the virial tensor.
314 * We could also determine the virial contribution directly,
315 * which would be cheaper here, but that would require extra
316 * communication for f_novirsum for with virtual sites
319 xf_z
[XX
] -= x
[i
][XX
]*F
;
320 xf_z
[YY
] -= x
[i
][YY
]*F
;
321 xf_z
[ZZ
] -= wall_z
[w
]*F
;
327 dvdlambda
+= (lam
== 0 ? -1 : 1)*Vlambda
;
330 inc_nrnb(nrnb
, eNR_WALLS
, md
->homenr
);
333 for (int i
= 0; i
< DIM
; i
++)
335 fr
->vir_wall_z
[i
] = -0.5*xf_z
[i
];