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-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
39 #include "vsite_parm.h"
47 #include "gromacs/gmxpreprocess/add_par.h"
48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
49 #include "gromacs/gmxpreprocess/grompp_impl.h"
50 #include "gromacs/gmxpreprocess/notset.h"
51 #include "gromacs/gmxpreprocess/toputil.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/mdtypes/md_enums.h"
56 #include "gromacs/topology/ifunc.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/utility/arrayref.h"
59 #include "gromacs/utility/cstringutil.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/gmxassert.h"
62 #include "gromacs/utility/smalloc.h"
63 #include "gromacs/utility/strconvert.h"
70 t_iatom
&ai() { return a
[0]; }
71 t_iatom
&aj() { return a
[1]; }
72 t_iatom
&ak() { return a
[2]; }
73 t_iatom
&al() { return a
[3]; }
76 struct VsiteBondParameter
78 VsiteBondParameter(int ftype
, const InteractionType
&type
)
79 : ftype_(ftype
), type_(type
)
82 const InteractionType
&type_
;
87 //! Function type for conversion.
89 //! The vsite parameters in a list.
90 std::vector
<VsiteBondParameter
> vSiteBondedParameters
;
93 using Atom2VsiteConnection
= std::vector
<int>;
95 static int vsite_bond_nrcheck(int ftype
)
99 if ((interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
| IF_ATYPE
)) || (ftype
== F_IDIHS
))
101 nrcheck
= NRAL(ftype
);
111 static void enter_bonded(int nratoms
, int *nrbonded
, t_mybonded
**bondeds
,
112 const InteractionType
&type
)
114 srenew(*bondeds
, *nrbonded
+1);
116 /* copy atom numbers */
117 gmx::ArrayRef
<const int> atoms
= type
.atoms();
118 GMX_RELEASE_ASSERT(nratoms
== atoms
.ssize(), "Size of atom array must much");
119 for (int j
= 0; j
< nratoms
; j
++)
121 (*bondeds
)[*nrbonded
].a
[j
] = atoms
[j
];
124 (*bondeds
)[*nrbonded
].c
= type
.c0();
129 static void get_bondeds(int nrat
, gmx::ArrayRef
<const int> atoms
,
130 gmx::ArrayRef
<const Atom2VsiteBond
> at2vb
,
131 int *nrbond
, t_mybonded
**bonds
,
132 int *nrang
, t_mybonded
**angles
,
133 int *nridih
, t_mybonded
**idihs
)
135 for (int k
= 0; k
< nrat
; k
++)
137 for (auto &vsite
: at2vb
[atoms
[k
]].vSiteBondedParameters
)
139 int ftype
= vsite
.ftype_
;
140 const InteractionType
&type
= vsite
.type_
;
141 int nrcheck
= vsite_bond_nrcheck(ftype
);
142 /* abuse nrcheck to see if we're adding bond, angle or idih */
145 case 2: enter_bonded(nrcheck
, nrbond
, bonds
, type
); break;
146 case 3: enter_bonded(nrcheck
, nrang
, angles
, type
); break;
147 case 4: enter_bonded(nrcheck
, nridih
, idihs
, type
); break;
153 static std::vector
<Atom2VsiteBond
>
154 make_at2vsitebond(int natoms
, gmx::ArrayRef
<InteractionTypeParameters
> plist
)
158 std::vector
<Atom2VsiteBond
> at2vb(natoms
);
161 for (int ftype
= 0; (ftype
< F_NRE
); ftype
++)
163 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
165 for (int i
= 0; (i
< gmx::ssize(plist
[ftype
])); i
++)
167 gmx::ArrayRef
<const int> atoms
= plist
[ftype
].interactionTypes
[i
].atoms();
168 for (int j
= 0; j
< NRAL(ftype
); j
++)
170 bVSI
[atoms
[j
]] = TRUE
;
176 for (int ftype
= 0; (ftype
< F_NRE
); ftype
++)
178 int nrcheck
= vsite_bond_nrcheck(ftype
);
181 for (int i
= 0; (i
< gmx::ssize(plist
[ftype
])); i
++)
183 gmx::ArrayRef
<const int> aa
= plist
[ftype
].interactionTypes
[i
].atoms();
184 for (int j
= 0; j
< nrcheck
; j
++)
188 at2vb
[aa
[j
]].vSiteBondedParameters
.emplace_back(ftype
, plist
[ftype
].interactionTypes
[i
]);
199 static std::vector
<Atom2VsiteConnection
>
200 make_at2vsitecon(int natoms
, gmx::ArrayRef
<InteractionTypeParameters
> plist
)
202 std::vector
<bool> bVSI(natoms
);
203 std::vector
<Atom2VsiteConnection
> at2vc(natoms
);
205 for (int ftype
= 0; (ftype
< F_NRE
); ftype
++)
207 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
209 for (int i
= 0; (i
< gmx::ssize(plist
[ftype
])); i
++)
211 gmx::ArrayRef
<const int> atoms
= plist
[ftype
].interactionTypes
[i
].atoms();
212 for (int j
= 0; j
< NRAL(ftype
); j
++)
214 bVSI
[atoms
[j
]] = TRUE
;
220 for (int ftype
= 0; (ftype
< F_NRE
); ftype
++)
222 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
224 for (int i
= 0; (i
< gmx::ssize(plist
[ftype
])); i
++)
226 int ai
= plist
[ftype
].interactionTypes
[i
].ai();
227 int aj
= plist
[ftype
].interactionTypes
[i
].aj();
228 if (bVSI
[ai
] && bVSI
[aj
])
230 /* Store forward direction */
231 at2vc
[ai
].emplace_back(aj
);
232 /* Store backward direction */
233 at2vc
[aj
].emplace_back(ai
);
242 static void print_bad(FILE *fp
,
243 int nrbond
, t_mybonded
*bonds
,
244 int nrang
, t_mybonded
*angles
,
245 int nridih
, t_mybonded
*idihs
)
249 fprintf(fp
, "bonds:");
250 for (int i
= 0; i
< nrbond
; i
++)
252 fprintf(fp
, " %d-%d (%g)",
253 bonds
[i
].ai()+1, bonds
[i
].aj()+1, bonds
[i
].c
);
259 fprintf(fp
, "angles:");
260 for (int i
= 0; i
< nrang
; i
++)
262 fprintf(fp
, " %d-%d-%d (%g)",
263 angles
[i
].ai()+1, angles
[i
].aj()+1,
264 angles
[i
].ak()+1, angles
[i
].c
);
270 fprintf(fp
, "idihs:");
271 for (int i
= 0; i
< nridih
; i
++)
273 fprintf(fp
, " %d-%d-%d-%d (%g)",
274 idihs
[i
].ai()+1, idihs
[i
].aj()+1,
275 idihs
[i
].ak()+1, idihs
[i
].al()+1, idihs
[i
].c
);
281 static void printInteractionType(FILE *fp
, int ftype
, int i
, const InteractionType
&type
)
284 static int prev_ftype
= NOTSET
;
285 static int prev_i
= NOTSET
;
287 if ( (ftype
!= prev_ftype
) || (i
!= prev_i
) )
293 fprintf(fp
, "(%d) plist[%s].param[%d]",
294 pass
, interaction_function
[ftype
].name
, i
);
295 gmx::ArrayRef
<const real
> forceParam
= type
.forceParam();
296 for (int j
= 0; j
< NRFP(ftype
); j
++)
298 fprintf(fp
, ".c[%d]=%g ", j
, forceParam
[j
]);
304 static real
get_bond_length(int nrbond
, t_mybonded bonds
[],
305 t_iatom ai
, t_iatom aj
)
311 for (i
= 0; i
< nrbond
&& (bondlen
== NOTSET
); i
++)
313 /* check both ways */
314 if ( ( (ai
== bonds
[i
].ai()) && (aj
== bonds
[i
].aj()) ) ||
315 ( (ai
== bonds
[i
].aj()) && (aj
== bonds
[i
].ai()) ) )
317 bondlen
= bonds
[i
].c
; /* note: bonds[i].c might be NOTSET */
323 static real
get_angle(int nrang
, t_mybonded angles
[],
324 t_iatom ai
, t_iatom aj
, t_iatom ak
)
330 for (i
= 0; i
< nrang
&& (angle
== NOTSET
); i
++)
332 /* check both ways */
333 if ( ( (ai
== angles
[i
].ai()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ak()) ) ||
334 ( (ai
== angles
[i
].ak()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ai()) ) )
336 angle
= DEG2RAD
*angles
[i
].c
;
342 static const char *get_atomtype_name_AB(t_atom
*atom
, PreprocessingAtomTypes
*atypes
)
344 const char* name
= atypes
->atomNameFromAtomType(atom
->type
);
346 /* When using the decoupling option, atom types are changed
347 * to decoupled for the non-bonded interactions, but the virtual
348 * sites constructions should be based on the "normal" interactions.
349 * So we return the state B atom type name if the state A atom
350 * type is the decoupled one. We should actually check for the atom
351 * type number, but that's not passed here. So we check for
352 * the decoupled atom type name. This should not cause trouble
353 * as this code is only used for topologies with v-sites without
354 * parameters generated by pdb2gmx.
356 if (strcmp(name
, "decoupled") == 0)
358 name
= atypes
->atomNameFromAtomType(atom
->typeB
);
364 static bool calc_vsite3_param(PreprocessingAtomTypes
*atypes
,
365 InteractionType
*param
, t_atoms
*at
,
366 int nrbond
, t_mybonded
*bonds
,
367 int nrang
, t_mybonded
*angles
)
369 /* i = virtual site | ,k
370 * j = 1st bonded heavy atom | i-j
371 * k,l = 2nd bonded atoms | `l
375 real bjk
, bjl
, a
= -1, b
= -1;
376 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
377 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
379 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atypes
), "MNH", 3) == 0) &&
380 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atypes
), "MNH", 3) == 0) ) ||
381 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atypes
), "MCH3", 4) == 0) &&
382 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atypes
), "MCH3", 4) == 0) );
384 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
385 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
386 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
389 /* now we get some XH2/XH3 group specific construction */
390 /* note: we call the heavy atom 'C' and the X atom 'N' */
391 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, dM
, rM
;
394 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
395 bError
= bError
|| (bjk
!= bjl
);
397 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
398 aN
= std::max(param
->ak(), param
->al())+1;
400 /* get common bonds */
401 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
403 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
404 bError
= bError
|| (bMM
== NOTSET
) || (bCN
== NOTSET
);
406 /* calculate common things */
408 dM
= std::sqrt( gmx::square(bCM
) - gmx::square(rM
) );
410 /* are we dealing with the X atom? */
411 if (param
->ai() == aN
)
413 /* this is trivial */
414 a
= b
= 0.5 * bCN
/dM
;
419 /* get other bondlengths and angles: */
420 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
421 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
422 bError
= bError
|| (bNH
== NOTSET
) || (aCNH
== NOTSET
);
425 dH
= bCN
- bNH
* std::cos(aCNH
);
426 rH
= bNH
* std::sin(aCNH
);
428 a
= 0.5 * ( dH
/dM
+ rH
/rM
);
429 b
= 0.5 * ( dH
/dM
- rH
/rM
);
434 gmx_fatal(FARGS
, "calc_vsite3_param not implemented for the general case "
435 "(atom %d)", param
->ai()+1);
437 param
->setForceParameter(0, a
);
438 param
->setForceParameter(1, b
);
443 static bool calc_vsite3fd_param(InteractionType
*param
,
444 int nrbond
, t_mybonded
*bonds
,
445 int nrang
, t_mybonded
*angles
)
447 /* i = virtual site | ,k
448 * j = 1st bonded heavy atom | i-j
449 * k,l = 2nd bonded atoms | `l
453 real bij
, bjk
, bjl
, aijk
, aijl
, rk
, rl
;
455 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
456 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
457 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
458 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
459 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
460 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) ||
461 (aijk
== NOTSET
) || (aijl
== NOTSET
);
463 rk
= bjk
* std::sin(aijk
);
464 rl
= bjl
* std::sin(aijl
);
465 param
->setForceParameter(0, rk
/ (rk
+ rl
));
466 param
->setForceParameter(1, -bij
);
471 static bool calc_vsite3fad_param(InteractionType
*param
,
472 int nrbond
, t_mybonded
*bonds
,
473 int nrang
, t_mybonded
*angles
)
475 /* i = virtual site |
476 * j = 1st bonded heavy atom | i-j
477 * k = 2nd bonded heavy atom | `k-l
478 * l = 3d bonded heavy atom |
481 bool bSwapParity
, bError
;
484 bSwapParity
= ( param
->c1() == -1 );
486 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
487 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
488 bError
= (bij
== NOTSET
) || (aijk
== NOTSET
);
490 param
->setForceParameter(1, bij
);
491 param
->setForceParameter(0, RAD2DEG
*aijk
);
495 param
->setForceParameter(0, 360 - param
->c0());
501 static bool calc_vsite3out_param(PreprocessingAtomTypes
*atypes
,
502 InteractionType
*param
, t_atoms
*at
,
503 int nrbond
, t_mybonded
*bonds
,
504 int nrang
, t_mybonded
*angles
)
506 /* i = virtual site | ,k
507 * j = 1st bonded heavy atom | i-j
508 * k,l = 2nd bonded atoms | `l
509 * NOTE: i is out of the j-k-l plane!
512 bool bXH3
, bError
, bSwapParity
;
513 real bij
, bjk
, bjl
, aijk
, aijl
, akjl
, pijk
, pijl
, a
, b
, c
;
515 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
516 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
518 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atypes
), "MNH", 3) == 0) &&
519 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atypes
), "MNH", 3) == 0) ) ||
520 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atypes
), "MCH3", 4) == 0) &&
521 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atypes
), "MCH3", 4) == 0) );
523 /* check if construction parity must be swapped */
524 bSwapParity
= ( param
->c1() == -1 );
526 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
527 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
528 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
531 /* now we get some XH3 group specific construction */
532 /* note: we call the heavy atom 'C' and the X atom 'N' */
533 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, rHx
, rHy
, dM
, rM
;
536 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
537 bError
= bError
|| (bjk
!= bjl
);
539 /* the X atom (C or N) in the XH3 group is the first after the masses: */
540 aN
= std::max(param
->ak(), param
->al())+1;
542 /* get all bondlengths and angles: */
543 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
545 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
546 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
547 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
549 (bMM
== NOTSET
) || (bCN
== NOTSET
) || (bNH
== NOTSET
) || (aCNH
== NOTSET
);
552 dH
= bCN
- bNH
* std::cos(aCNH
);
553 rH
= bNH
* std::sin(aCNH
);
554 /* we assume the H's are symmetrically distributed */
555 rHx
= rH
*std::cos(DEG2RAD
*30);
556 rHy
= rH
*std::sin(DEG2RAD
*30);
558 dM
= std::sqrt( gmx::square(bCM
) - gmx::square(rM
) );
559 a
= 0.5*( (dH
/dM
) - (rHy
/rM
) );
560 b
= 0.5*( (dH
/dM
) + (rHy
/rM
) );
566 /* this is the general construction */
568 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
569 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
570 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
571 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
573 (bij
== NOTSET
) || (aijk
== NOTSET
) || (aijl
== NOTSET
) || (akjl
== NOTSET
);
575 pijk
= std::cos(aijk
)*bij
;
576 pijl
= std::cos(aijl
)*bij
;
577 a
= ( pijk
+ (pijk
*std::cos(akjl
)-pijl
) * std::cos(akjl
) / gmx::square(std::sin(akjl
)) ) / bjk
;
578 b
= ( pijl
+ (pijl
*std::cos(akjl
)-pijk
) * std::cos(akjl
) / gmx::square(std::sin(akjl
)) ) / bjl
;
579 c
= -std::sqrt( gmx::square(bij
) -
580 ( gmx::square(pijk
) - 2*pijk
*pijl
*std::cos(akjl
) + gmx::square(pijl
) )
581 / gmx::square(std::sin(akjl
)) )
582 / ( bjk
*bjl
*std::sin(akjl
) );
585 param
->setForceParameter(0, a
);
586 param
->setForceParameter(1, b
);
589 param
->setForceParameter(2, -c
);
593 param
->setForceParameter(2, c
);
598 static bool calc_vsite4fd_param(InteractionType
*param
,
599 int nrbond
, t_mybonded
*bonds
,
600 int nrang
, t_mybonded
*angles
)
602 /* i = virtual site | ,k
603 * j = 1st bonded heavy atom | i-j-m
604 * k,l,m = 2nd bonded atoms | `l
608 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
, akjm
, akjl
;
609 real pk
, pl
, pm
, cosakl
, cosakm
, sinakl
, sinakm
, cl
, cm
;
611 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
612 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
613 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
614 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
615 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
616 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
617 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
618 akjm
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->am());
619 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
620 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
621 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
) || (akjm
== NOTSET
) ||
626 pk
= bjk
*std::sin(aijk
);
627 pl
= bjl
*std::sin(aijl
);
628 pm
= bjm
*std::sin(aijm
);
629 cosakl
= (std::cos(akjl
) - std::cos(aijk
)*std::cos(aijl
)) / (std::sin(aijk
)*std::sin(aijl
));
630 cosakm
= (std::cos(akjm
) - std::cos(aijk
)*std::cos(aijm
)) / (std::sin(aijk
)*std::sin(aijm
));
631 if (cosakl
< -1 || cosakl
> 1 || cosakm
< -1 || cosakm
> 1)
633 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
634 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
635 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fd for atom %d: "
636 "cosakl=%f, cosakm=%f\n", param
->ai()+1, cosakl
, cosakm
);
638 sinakl
= std::sqrt(1-gmx::square(cosakl
));
639 sinakm
= std::sqrt(1-gmx::square(cosakm
));
641 /* note: there is a '+' because of the way the sines are calculated */
642 cl
= -pk
/ ( pl
*cosakl
- pk
+ pl
*sinakl
*(pm
*cosakm
-pk
)/(pm
*sinakm
) );
643 cm
= -pk
/ ( pm
*cosakm
- pk
+ pm
*sinakm
*(pl
*cosakl
-pk
)/(pl
*sinakl
) );
645 param
->setForceParameter(0, cl
);
646 param
->setForceParameter(1, cm
);
647 param
->setForceParameter(2, -bij
);
655 calc_vsite4fdn_param(InteractionType
*param
,
656 int nrbond
, t_mybonded
*bonds
,
657 int nrang
, t_mybonded
*angles
)
659 /* i = virtual site | ,k
660 * j = 1st bonded heavy atom | i-j-m
661 * k,l,m = 2nd bonded atoms | `l
665 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
;
666 real pk
, pl
, pm
, a
, b
;
668 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
669 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
670 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
671 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
672 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
673 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
674 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
676 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
677 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
);
682 /* Calculate component of bond j-k along the direction i-j */
683 pk
= -bjk
*std::cos(aijk
);
685 /* Calculate component of bond j-l along the direction i-j */
686 pl
= -bjl
*std::cos(aijl
);
688 /* Calculate component of bond j-m along the direction i-j */
689 pm
= -bjm
*std::cos(aijm
);
691 if (fabs(pl
) < 1000*GMX_REAL_MIN
|| fabs(pm
) < 1000*GMX_REAL_MIN
)
693 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
694 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
695 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fdn for atom %d: "
696 "pl=%f, pm=%f\n", param
->ai()+1, pl
, pm
);
702 param
->setForceParameter(0, a
);
703 param
->setForceParameter(1, b
);
704 param
->setForceParameter(2, bij
);
712 int set_vsites(bool bVerbose
, t_atoms
*atoms
, PreprocessingAtomTypes
*atypes
,
713 gmx::ArrayRef
<InteractionTypeParameters
> plist
)
716 int nvsite
, nrbond
, nrang
, nridih
, nrset
;
725 /* Make a reverse list to avoid ninteractions^2 operations */
726 std::vector
<Atom2VsiteBond
> at2vb
= make_at2vsitebond(atoms
->nr
, plist
);
728 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
730 if (interaction_function
[ftype
].flags
& IF_VSITE
)
732 nvsite
+= plist
[ftype
].size();
734 if (ftype
== F_VSITEN
)
736 /* We don't calculate parameters for VSITEN */
742 for (auto ¶m
: plist
[ftype
].interactionTypes
)
744 /* check if all parameters are set */
746 gmx::ArrayRef
<const real
> forceParam
= param
.forceParam();
747 for (int j
= 0; (j
< NRFP(ftype
)) && bSet
; j
++)
749 bSet
= forceParam
[j
] != NOTSET
;
754 fprintf(debug
, "bSet=%s ", gmx::boolToString(bSet
));
755 printInteractionType(debug
, ftype
, i
, plist
[ftype
].interactionTypes
[i
]);
759 if (bVerbose
&& bFirst
)
761 fprintf(stderr
, "Calculating parameters for virtual sites\n");
765 nrbond
= nrang
= nridih
= 0;
770 /* now set the vsite parameters: */
771 get_bondeds(NRAL(ftype
), param
.atoms(), at2vb
,
772 &nrbond
, &bonds
, &nrang
, &angles
, &nridih
, &idihs
);
775 fprintf(debug
, "Found %d bonds, %d angles and %d idihs "
776 "for virtual site %d (%s)\n", nrbond
, nrang
, nridih
,
778 interaction_function
[ftype
].longname
);
779 print_bad(debug
, nrbond
, bonds
, nrang
, angles
, nridih
, idihs
);
785 calc_vsite3_param(atypes
, ¶m
, atoms
,
786 nrbond
, bonds
, nrang
, angles
);
790 calc_vsite3fd_param(¶m
,
791 nrbond
, bonds
, nrang
, angles
);
795 calc_vsite3fad_param(¶m
,
796 nrbond
, bonds
, nrang
, angles
);
800 calc_vsite3out_param(atypes
, ¶m
, atoms
,
801 nrbond
, bonds
, nrang
, angles
);
805 calc_vsite4fd_param(¶m
,
806 nrbond
, bonds
, nrang
, angles
);
810 calc_vsite4fdn_param(¶m
,
811 nrbond
, bonds
, nrang
, angles
);
814 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
816 interaction_function
[ftype
].longname
,
822 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
823 "for %s atom %d for this bonding configuration",
824 interaction_function
[ftype
].longname
,
839 void set_vsites_ptype(bool bVerbose
, gmx_moltype_t
*molt
)
845 fprintf(stderr
, "Setting particle type to V for virtual sites\n");
847 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
849 InteractionList
*il
= &molt
->ilist
[ftype
];
850 if (interaction_function
[ftype
].flags
& IF_VSITE
)
852 const int nra
= interaction_function
[ftype
].nratoms
;
853 const int nrd
= il
->size();
854 gmx::ArrayRef
<const int> ia
= il
->iatoms
;
858 fprintf(stderr
, "doing %d %s virtual sites\n",
859 (nrd
/ (nra
+1)), interaction_function
[ftype
].longname
);
862 for (i
= 0; (i
< nrd
); )
864 /* The virtual site */
865 int avsite
= ia
[i
+ 1];
866 molt
->atoms
.atom
[avsite
].ptype
= eptVSite
;
879 static void check_vsite_constraints(gmx::ArrayRef
<InteractionTypeParameters
> plist
,
880 int cftype
, const int vsite_type
[])
883 for (const auto ¶m
: plist
[cftype
].interactionTypes
)
885 gmx::ArrayRef
<const int> atoms
= param
.atoms();
886 for (int k
= 0; k
< 2; k
++)
889 if (vsite_type
[atom
] != NOTSET
)
891 fprintf(stderr
, "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)\n",
892 param
.ai()+1, param
.aj()+1, atom
+1);
899 gmx_fatal(FARGS
, "There were %d virtual sites involved in constraints", n
);
903 static void clean_vsite_bonds(gmx::ArrayRef
<InteractionTypeParameters
> plist
, t_pindex pindex
[],
904 int cftype
, const int vsite_type
[])
907 int nconverted
, nremoved
;
909 bool bKeep
, bRemove
, bAllFD
;
910 InteractionTypeParameters
*ps
;
912 if (cftype
== F_CONNBONDS
)
917 ps
= &(plist
[cftype
]);
921 for (auto bond
= ps
->interactionTypes
.begin(); bond
!= ps
->interactionTypes
.end(); )
924 const int *first_atoms
= nullptr;
929 /* check if all virtual sites are constructed from the same atoms */
931 gmx::ArrayRef
<const int> atoms
= bond
->atoms();
932 for (int k
= 0; (k
< 2) && !bKeep
&& !bRemove
; k
++)
934 /* for all atoms in the bond */
936 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
939 bool bThisFD
= ( (pindex
[atom
].ftype
== F_VSITE3FD
) ||
940 (pindex
[atom
].ftype
== F_VSITE3FAD
) ||
941 (pindex
[atom
].ftype
== F_VSITE4FD
) ||
942 (pindex
[atom
].ftype
== F_VSITE4FDN
) );
943 bool bThisOUT
= ( (pindex
[atom
].ftype
== F_VSITE3OUT
) &&
944 ((interaction_function
[cftype
].flags
& IF_CONSTRAINT
) != 0u) );
945 bAllFD
= bAllFD
&& bThisFD
;
946 if (bThisFD
|| bThisOUT
)
948 oatom
= atoms
[1-k
]; /* the other atom */
949 if (vsite_type
[oatom
] == NOTSET
&&
950 oatom
== plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].aj())
952 /* if the other atom isn't a vsite, and it is AI */
962 /* TODO This fragment, and corresponding logic in
963 clean_vsite_angles and clean_vsite_dihs should
964 be refactored into a common function */
967 /* if this is the first vsite we encounter then
968 store construction atoms */
969 /* TODO This would be nicer to implement with
970 a C++ "vector view" class" with an
971 STL-container-like interface. */
972 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
973 first_atoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
977 GMX_ASSERT(vsnral
!= 0, "nvsite > 1 must have vsnral != 0");
978 GMX_ASSERT(first_atoms
!= nullptr, "nvsite > 1 must have first_atoms != NULL");
979 /* if it is not the first then
980 check if this vsite is constructed from the same atoms */
981 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
983 for (int m
= 0; (m
< vsnral
) && !bKeep
; m
++)
987 bool bPresent
= false;
988 atoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
989 for (int n
= 0; (n
< vsnral
) && !bPresent
; n
++)
991 if (atoms
[m
] == first_atoms
[n
])
1017 /* if we have no virtual sites in this bond, keep it */
1023 /* TODO This loop and the corresponding loop in
1024 check_vsite_angles should be refactored into a common
1026 /* check if all non-vsite atoms are used in construction: */
1027 bool bFirstTwo
= true;
1028 for (int k
= 0; (k
< 2) && !bKeep
; k
++) /* for all atoms in the bond */
1030 int atom
= atoms
[k
];
1031 if (vsite_type
[atom
] == NOTSET
)
1034 for (int m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1036 GMX_ASSERT(first_atoms
!= nullptr, "If we've seen a vsite before, we know what its first atom index was");
1038 if (atom
== first_atoms
[m
])
1041 bFirstTwo
= bFirstTwo
&& m
< 2;
1051 if (!( bAllFD
&& bFirstTwo
) )
1053 /* Two atom bonded interactions include constraints.
1054 * We need to remove constraints between vsite pairs that have
1055 * a fixed distance due to being constructed from the same
1056 * atoms, since this can be numerically unstable.
1058 for (int m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1060 at1
= first_atoms
[m
];
1061 at2
= first_atoms
[(m
+1) % vsnral
];
1062 bool bPresent
= false;
1063 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1065 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1067 for (auto entry
= plist
[ftype
].interactionTypes
.begin(); (entry
!= plist
[ftype
].interactionTypes
.end()) && !bPresent
; entry
++)
1069 /* all constraints until one matches */
1070 bPresent
= ( ( (entry
->ai() == at1
) &&
1071 (entry
->aj() == at2
) ) ||
1072 ( (entry
->ai() == at2
) &&
1073 (entry
->aj() == at1
) ) );
1089 else if (IS_CHEMBOND(cftype
))
1091 plist
[F_CONNBONDS
].interactionTypes
.emplace_back(*bond
);
1092 bond
= ps
->interactionTypes
.erase(bond
);
1097 bond
= ps
->interactionTypes
.erase(bond
);
1104 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %zu left\n",
1105 nremoved
, interaction_function
[cftype
].longname
, ps
->size());
1109 fprintf(stderr
, "Converted %4d %15ss with virtual sites to connections, %zu left\n",
1110 nconverted
, interaction_function
[cftype
].longname
, ps
->size());
1114 fprintf(stderr
, "Warning: removed %d %ss with vsite with %s construction\n"
1115 " This vsite construction does not guarantee constant "
1117 " If the constructions were generated by pdb2gmx ignore "
1119 nOut
, interaction_function
[cftype
].longname
,
1120 interaction_function
[F_VSITE3OUT
].longname
);
1124 static void clean_vsite_angles(gmx::ArrayRef
<InteractionTypeParameters
> plist
, t_pindex pindex
[],
1125 int cftype
, const int vsite_type
[],
1126 gmx::ArrayRef
<const Atom2VsiteConnection
> at2vc
)
1129 InteractionTypeParameters
*ps
;
1131 ps
= &(plist
[cftype
]);
1132 int oldSize
= ps
->size();
1133 for (auto angle
= ps
->interactionTypes
.begin(); angle
!= ps
->interactionTypes
.end(); )
1136 const int *first_atoms
= nullptr;
1139 bool bAll3FAD
= true;
1140 /* check if all virtual sites are constructed from the same atoms */
1142 gmx::ArrayRef
<const int> atoms
= angle
->atoms();
1143 for (int k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1145 int atom
= atoms
[k
];
1146 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1149 bAll3FAD
= bAll3FAD
&& (pindex
[atom
].ftype
== F_VSITE3FAD
);
1152 /* store construction atoms of first vsite */
1153 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1154 first_atoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
1158 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1159 GMX_ASSERT(first_atoms
!= nullptr, "If we've seen a vsite before, we know what its first atom index was");
1160 /* check if this vsite is constructed from the same atoms */
1161 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1163 for (int m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1165 const int *subAtoms
;
1167 bool bPresent
= false;
1168 subAtoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
1169 for (int n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1171 if (subAtoms
[m
] == first_atoms
[n
])
1190 /* keep all angles with no virtual sites in them or
1191 with virtual sites with more than 3 constr. atoms */
1192 if (nvsite
== 0 && vsnral
> 3)
1197 /* check if all non-vsite atoms are used in construction: */
1198 bool bFirstTwo
= true;
1199 for (int k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1202 if (vsite_type
[atom
] == NOTSET
)
1205 for (int m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1207 GMX_ASSERT(first_atoms
!= nullptr, "If we've seen a vsite before, we know what its first atom index was");
1209 if (atom
== first_atoms
[m
])
1212 bFirstTwo
= bFirstTwo
&& m
< 2;
1222 if (!( bAll3FAD
&& bFirstTwo
) )
1224 /* check if all constructing atoms are constrained together */
1225 for (int m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1227 at1
= first_atoms
[m
];
1228 at2
= first_atoms
[(m
+1) % vsnral
];
1229 bool bPresent
= false;
1230 auto found
= std::find(at2vc
[at1
].begin(), at2vc
[at1
].end(), at2
);
1231 if (found
!= at2vc
[at1
].end())
1248 angle
= ps
->interactionTypes
.erase(angle
);
1252 if (oldSize
!= gmx::ssize(*ps
))
1254 fprintf(stderr
, "Removed %4zu %15ss with virtual sites, %zu left\n",
1255 oldSize
-ps
->size(), interaction_function
[cftype
].longname
, ps
->size());
1259 static void clean_vsite_dihs(gmx::ArrayRef
<InteractionTypeParameters
> plist
, t_pindex pindex
[],
1260 int cftype
, const int vsite_type
[])
1262 InteractionTypeParameters
*ps
;
1264 ps
= &(plist
[cftype
]);
1266 int oldSize
= ps
->size();
1267 for (auto dih
= ps
->interactionTypes
.begin(); dih
!= ps
->interactionTypes
.end(); )
1270 const int *first_atoms
= nullptr;
1273 gmx::ArrayRef
<const int> atoms
= dih
->atoms();
1275 /* check if all virtual sites are constructed from the same atoms */
1277 for (int k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1280 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1284 /* store construction atoms of first vsite */
1285 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1286 first_atoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
1290 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1291 GMX_ASSERT(first_atoms
!= nullptr, "If we've seen a vsite before, we know what its first atom index was");
1292 /* check if this vsite is constructed from the same atoms */
1293 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1295 for (int m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1297 const int *subAtoms
;
1299 bool bPresent
= false;
1300 subAtoms
= plist
[pindex
[atom
].ftype
].interactionTypes
[pindex
[atom
].parnr
].atoms().data() + 1;
1301 for (int n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1303 if (subAtoms
[m
] == first_atoms
[n
])
1315 /* TODO clean_site_bonds and _angles do this increment
1316 at the top of the loop. Refactor this for
1322 /* keep all dihedrals with no virtual sites in them */
1328 /* check if all atoms in dihedral are either virtual sites, or used in
1329 construction of virtual sites. If so, keep it, if not throw away: */
1330 for (int k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1332 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1333 GMX_ASSERT(first_atoms
!= nullptr, "If we've seen a vsite before, we know what its first atom index was");
1335 if (vsite_type
[atom
] == NOTSET
)
1337 /* vsnral will be set here, we don't get here with nvsite==0 */
1339 for (int m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1341 if (atom
== first_atoms
[m
])
1359 dih
= ps
->interactionTypes
.erase(dih
);
1364 if (oldSize
!= gmx::ssize(*ps
))
1366 fprintf(stderr
, "Removed %4zu %15ss with virtual sites, %zu left\n",
1367 oldSize
-ps
->size(), interaction_function
[cftype
].longname
, ps
->size());
1371 void clean_vsite_bondeds(gmx::ArrayRef
<InteractionTypeParameters
> plist
, int natoms
, bool bRmVSiteBds
)
1376 std::vector
<Atom2VsiteConnection
> at2vc
;
1378 pindex
= nullptr; /* avoid warnings */
1379 /* make vsite_type array */
1380 snew(vsite_type
, natoms
);
1381 for (int i
= 0; i
< natoms
; i
++)
1383 vsite_type
[i
] = NOTSET
;
1386 for (int ftype
= 0; ftype
< F_NRE
; ftype
++)
1388 if (interaction_function
[ftype
].flags
& IF_VSITE
)
1390 nvsite
+= plist
[ftype
].size();
1392 while (i
< gmx::ssize(plist
[ftype
]))
1394 vsite
= plist
[ftype
].interactionTypes
[i
].ai();
1395 if (vsite_type
[vsite
] == NOTSET
)
1397 vsite_type
[vsite
] = ftype
;
1401 gmx_fatal(FARGS
, "multiple vsite constructions for atom %d", vsite
+1);
1403 if (ftype
== F_VSITEN
)
1405 while (i
< gmx::ssize(plist
[ftype
]) && plist
[ftype
].interactionTypes
[i
].ai() == vsite
)
1418 /* the rest only if we have virtual sites: */
1421 fprintf(stderr
, "Cleaning up constraints %swith virtual sites\n",
1422 bRmVSiteBds
? "and constant bonded interactions " : "");
1424 /* Make a reverse list to avoid ninteractions^2 operations */
1425 at2vc
= make_at2vsitecon(natoms
, plist
);
1427 snew(pindex
, natoms
);
1428 for (int ftype
= 0; ftype
< F_NRE
; ftype
++)
1430 /* Here we skip VSITEN. In neary all practical use cases this
1431 * is not an issue, since VSITEN is intended for constructing
1432 * additional interaction sites, not for replacing normal atoms
1433 * with bonded interactions. Thus we do not expect constant
1434 * bonded interactions. If VSITEN does get used with constant
1435 * bonded interactions, these are not removed which only leads
1436 * to very minor extra computation and constant energy.
1437 * The only problematic case is onstraints between VSITEN
1438 * constructions with fixed distance (which is anyhow useless).
1439 * This will generate a fatal error in check_vsite_constraints.
1441 if ((interaction_function
[ftype
].flags
& IF_VSITE
) &&
1444 for (int parnr
= 0; (parnr
< gmx::ssize(plist
[ftype
])); parnr
++)
1446 int k
= plist
[ftype
].interactionTypes
[parnr
].ai();
1447 pindex
[k
].ftype
= ftype
;
1448 pindex
[k
].parnr
= parnr
;
1453 /* remove interactions that include virtual sites */
1454 for (int ftype
= 0; ftype
< F_NRE
; ftype
++)
1456 if ( ( ( interaction_function
[ftype
].flags
& IF_BOND
) && bRmVSiteBds
) ||
1457 ( interaction_function
[ftype
].flags
& IF_CONSTRAINT
) )
1459 if (interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
) )
1461 clean_vsite_bonds (plist
, pindex
, ftype
, vsite_type
);
1463 else if (interaction_function
[ftype
].flags
& IF_ATYPE
)
1465 clean_vsite_angles(plist
, pindex
, ftype
, vsite_type
, at2vc
);
1467 else if ( (ftype
== F_PDIHS
) || (ftype
== F_IDIHS
) )
1469 clean_vsite_dihs (plist
, pindex
, ftype
, vsite_type
);
1473 /* check that no remaining constraints include virtual sites */
1474 for (int ftype
= 0; ftype
< F_NRE
; ftype
++)
1476 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1478 check_vsite_constraints(plist
, ftype
, vsite_type
);