#include "gromacs/gmxpreprocess/vsite_parm.h"
#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
-#include "gromacs/mdlib/genborn.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/pbcutil/pbc.h"
}
-static int
-find_gb_bondlength(t_params *plist, int ai, int aj, real *length)
-{
- int i, j, a1, a2;
-
- int found = 0;
- int status;
-
- for (i = 0; i < F_NRE && !found; i++)
- {
- if (IS_CHEMBOND(i))
- {
- for (j = 0; j < plist[i].nr; j++)
- {
- a1 = plist[i].param[j].a[0];
- a2 = plist[i].param[j].a[1];
-
- if ( (a1 == ai && a2 == aj) || (a1 == aj && a2 == ai))
- {
- /* Equilibrium bond distance */
- *length = plist[i].param[j].c[0];
- found = 1;
- }
- }
- }
- }
- status = !found;
-
- return status;
-}
-
-
-static int
-find_gb_anglelength(t_params *plist, int ai, int ak, real *length)
-{
- int i, j, a1, a2, a3;
- real r12, r23, a123;
- int found = 0;
- int status, status1, status2;
-
- r12 = r23 = 0;
-
- for (i = 0; i < F_NRE && !found; i++)
- {
- if (IS_ANGLE(i))
- {
- for (j = 0; j < plist[i].nr; j++)
- {
- a1 = plist[i].param[j].a[0];
- a2 = plist[i].param[j].a[1];
- a3 = plist[i].param[j].a[2];
-
- /* We dont care what the middle atom is, but use it below */
- if ( (a1 == ai && a3 == ak) || (a1 == ak && a3 == ai) )
- {
- /* Equilibrium bond distance */
- a123 = plist[i].param[j].c[0];
- /* Use middle atom to find reference distances r12 and r23 */
- status1 = find_gb_bondlength(plist, a1, a2, &r12);
- status2 = find_gb_bondlength(plist, a2, a3, &r23);
-
- if (status1 == 0 && status2 == 0)
- {
- /* cosine theorem to get r13 */
- *length = std::sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG)));
- found = 1;
- }
- }
- }
- }
- }
- status = !found;
-
- return status;
-}
-
-static int
-generate_gb_exclusion_interactions(t_molinfo *mi, gpp_atomtype_t atype, t_nextnb *nnb)
-{
- int j, n, ai, aj, ti, tj;
- int ftype;
- t_param param;
- t_params * plist;
- t_atoms * at;
- real radiusi, radiusj;
- real gb_radiusi, gb_radiusj;
- real param_c2, param_c4;
- real distance;
-
- plist = mi->plist;
- at = &mi->atoms;
-
- for (n = 1; n <= nnb->nrex; n++)
- {
- switch (n)
- {
- case 1:
- ftype = F_GB12;
- param_c2 = STILL_P2;
- param_c4 = 0.8875;
- break;
- case 2:
- ftype = F_GB13;
- param_c2 = STILL_P3;
- param_c4 = 0.3516;
- break;
- default:
- /* Put all higher-order exclusions into 1,4 list so we dont miss them */
- ftype = F_GB14;
- param_c2 = STILL_P3;
- param_c4 = 0.3516;
- break;
- }
-
- for (ai = 0; ai < nnb->nr; ai++)
- {
- ti = at->atom[ai].type;
- radiusi = get_atomtype_radius(ti, atype);
- gb_radiusi = get_atomtype_gb_radius(ti, atype);
-
- for (j = 0; j < nnb->nrexcl[ai][n]; j++)
- {
- aj = nnb->a[ai][n][j];
-
- /* Only add the interactions once */
- if (aj > ai)
- {
- tj = at->atom[aj].type;
- radiusj = get_atomtype_radius(tj, atype);
- gb_radiusj = get_atomtype_gb_radius(tj, atype);
-
- /* There is an exclusion of type "ftype" between atoms ai and aj */
- param.a[0] = ai;
- param.a[1] = aj;
-
- /* Reference distance, not used for 1-4 interactions */
- switch (ftype)
- {
- case F_GB12:
- if (find_gb_bondlength(plist, ai, aj, &distance) != 0)
- {
- gmx_fatal(FARGS, "Cannot find bond length for atoms %d-%d", ai, aj);
- }
- break;
- case F_GB13:
- if (find_gb_anglelength(plist, ai, aj, &distance) != 0)
- {
- gmx_fatal(FARGS, "Cannot find length for atoms %d-%d involved in angle", ai, aj);
- }
- break;
- default:
- distance = -1;
- break;
- }
- /* Assign GB parameters */
- /* Sum of radii */
- param.c[0] = radiusi+radiusj;
- /* Reference distance distance */
- param.c[1] = distance;
- /* Still parameter */
- param.c[2] = param_c2;
- /* GB radius */
- param.c[3] = gb_radiusi+gb_radiusj;
- /* Parameter */
- param.c[4] = param_c4;
-
- /* Add it to the parameter list */
- add_param_to_list(&plist[ftype], ¶m);
- }
- }
- }
- }
- return 0;
-}
-
-
static void make_atoms_sys(int nmolb, const gmx_molblock_t *molb,
const t_molinfo *molinfo,
t_atoms *atoms)
int *nmolblock,
gmx_molblock_t **molblock,
gmx_bool bFEP,
- gmx_bool bGenborn,
gmx_bool bZero,
gmx_bool usingFullRangeElectrostatics,
warninp_t wi)
*/
case d_implicit_genborn_params:
- push_gb_params(atype, pline, wi);
+ // Skip this line, so old topologies with
+ // GB parameters can be read.
break;
case d_implicit_surface_params:
- gmx_fatal(FARGS, "Implicit surface directive not supported yet.");
+ // Skip this line, so that any topologies
+ // with surface parameters can be read
+ // (even though these were never formally
+ // supported).
break;
case d_cmaptypes:
- /* nnb contains information about first,2nd,3rd bonded neighbors.
- * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
- */
- if (bGenborn == TRUE)
- {
- generate_gb_exclusion_interactions(mi0, atype, &nnb);
- }
-
done_nnb(&nnb);
if (bCouple)
const t_inputrec *ir,
int *nmolblock,
gmx_molblock_t **molblock,
- gmx_bool bGenborn,
warninp_t wi)
{
/* Tmpfile might contain a long path */
nrmols, molinfo, intermolecular_interactions,
plist, combination_rule, repulsion_power,
opts, fudgeQQ, nmolblock, molblock,
- ir->efep != efepNO, bGenborn, bZero,
+ ir->efep != efepNO, bZero,
EEL_FULL(ir->coulombtype), wi);
if ((*combination_rule != eCOMB_GEOMETRIC) &&
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff