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, 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"
48 #include "gromacs/gmxpreprocess/add_par.h"
49 #include "gromacs/gmxpreprocess/notset.h"
50 #include "gromacs/gmxpreprocess/resall.h"
51 #include "gromacs/gmxpreprocess/toputil.h"
52 #include "gromacs/legacyheaders/names.h"
53 #include "gromacs/legacyheaders/types/ifunc.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/gmxassert.h"
59 #include "gromacs/utility/smalloc.h"
64 t_iatom
&ai() { return a
[0]; }
65 t_iatom
&aj() { return a
[1]; }
66 t_iatom
&ak() { return a
[2]; }
67 t_iatom
&al() { return a
[3]; }
78 vsitebondparam_t
*vsbp
;
86 static int vsite_bond_nrcheck(int ftype
)
90 if ((interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
| IF_ATYPE
)) || (ftype
== F_IDIHS
))
92 nrcheck
= NRAL(ftype
);
102 static void enter_bonded(int nratoms
, int *nrbonded
, t_mybonded
**bondeds
,
107 srenew(*bondeds
, *nrbonded
+1);
109 /* copy atom numbers */
110 for (j
= 0; j
< nratoms
; j
++)
112 (*bondeds
)[*nrbonded
].a
[j
] = param
->a
[j
];
115 (*bondeds
)[*nrbonded
].c
= param
->c0();
120 static void get_bondeds(int nrat
, t_iatom atoms
[],
121 at2vsitebond_t
*at2vb
,
122 int *nrbond
, t_mybonded
**bonds
,
123 int *nrang
, t_mybonded
**angles
,
124 int *nridih
, t_mybonded
**idihs
)
126 int k
, i
, ftype
, nrcheck
;
129 for (k
= 0; k
< nrat
; k
++)
131 for (i
= 0; i
< at2vb
[atoms
[k
]].nr
; i
++)
133 ftype
= at2vb
[atoms
[k
]].vsbp
[i
].ftype
;
134 param
= at2vb
[atoms
[k
]].vsbp
[i
].param
;
135 nrcheck
= vsite_bond_nrcheck(ftype
);
136 /* abuse nrcheck to see if we're adding bond, angle or idih */
139 case 2: enter_bonded(nrcheck
, nrbond
, bonds
, param
); break;
140 case 3: enter_bonded(nrcheck
, nrang
, angles
, param
); break;
141 case 4: enter_bonded(nrcheck
, nridih
, idihs
, param
); break;
147 static at2vsitebond_t
*make_at2vsitebond(int natoms
, t_params plist
[])
150 int ftype
, i
, j
, nrcheck
, nr
;
152 at2vsitebond_t
*at2vb
;
157 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
159 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
161 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
163 for (j
= 0; j
< NRAL(ftype
); j
++)
165 bVSI
[plist
[ftype
].param
[i
].a
[j
]] = TRUE
;
171 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
173 nrcheck
= vsite_bond_nrcheck(ftype
);
176 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
178 aa
= plist
[ftype
].param
[i
].a
;
179 for (j
= 0; j
< nrcheck
; j
++)
183 nr
= at2vb
[aa
[j
]].nr
;
186 srenew(at2vb
[aa
[j
]].vsbp
, nr
+10);
188 at2vb
[aa
[j
]].vsbp
[nr
].ftype
= ftype
;
189 at2vb
[aa
[j
]].vsbp
[nr
].param
= &plist
[ftype
].param
[i
];
201 static void done_at2vsitebond(int natoms
, at2vsitebond_t
*at2vb
)
205 for (i
= 0; i
< natoms
; i
++)
209 sfree(at2vb
[i
].vsbp
);
215 static at2vsitecon_t
*make_at2vsitecon(int natoms
, t_params plist
[])
218 int ftype
, i
, j
, ai
, aj
, nr
;
219 at2vsitecon_t
*at2vc
;
224 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
226 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
228 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
230 for (j
= 0; j
< NRAL(ftype
); j
++)
232 bVSI
[plist
[ftype
].param
[i
].a
[j
]] = TRUE
;
238 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
240 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
242 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
244 ai
= plist
[ftype
].param
[i
].ai();
245 aj
= plist
[ftype
].param
[i
].aj();
246 if (bVSI
[ai
] && bVSI
[aj
])
248 /* Store forward direction */
252 srenew(at2vc
[ai
].aj
, nr
+10);
254 at2vc
[ai
].aj
[nr
] = aj
;
256 /* Store backward direction */
260 srenew(at2vc
[aj
].aj
, nr
+10);
262 at2vc
[aj
].aj
[nr
] = ai
;
273 static void done_at2vsitecon(int natoms
, at2vsitecon_t
*at2vc
)
277 for (i
= 0; i
< natoms
; i
++)
288 static void print_bad(FILE *fp
,
289 int nrbond
, t_mybonded
*bonds
,
290 int nrang
, t_mybonded
*angles
,
291 int nridih
, t_mybonded
*idihs
)
297 fprintf(fp
, "bonds:");
298 for (i
= 0; i
< nrbond
; i
++)
300 fprintf(fp
, " %d-%d (%g)",
301 bonds
[i
].ai()+1, bonds
[i
].aj()+1, bonds
[i
].c
);
307 fprintf(fp
, "angles:");
308 for (i
= 0; i
< nrang
; i
++)
310 fprintf(fp
, " %d-%d-%d (%g)",
311 angles
[i
].ai()+1, angles
[i
].aj()+1,
312 angles
[i
].ak()+1, angles
[i
].c
);
318 fprintf(fp
, "idihs:");
319 for (i
= 0; i
< nridih
; i
++)
321 fprintf(fp
, " %d-%d-%d-%d (%g)",
322 idihs
[i
].ai()+1, idihs
[i
].aj()+1,
323 idihs
[i
].ak()+1, idihs
[i
].al()+1, idihs
[i
].c
);
329 static void print_param(FILE *fp
, int ftype
, int i
, t_param
*param
)
332 static int prev_ftype
= NOTSET
;
333 static int prev_i
= NOTSET
;
336 if ( (ftype
!= prev_ftype
) || (i
!= prev_i
) )
342 fprintf(fp
, "(%d) plist[%s].param[%d]",
343 pass
, interaction_function
[ftype
].name
, i
);
344 for (j
= 0; j
< NRFP(ftype
); j
++)
346 fprintf(fp
, ".c[%d]=%g ", j
, param
->c
[j
]);
352 static real
get_bond_length(int nrbond
, t_mybonded bonds
[],
353 t_iatom ai
, t_iatom aj
)
359 for (i
= 0; i
< nrbond
&& (bondlen
== NOTSET
); i
++)
361 /* check both ways */
362 if ( ( (ai
== bonds
[i
].ai()) && (aj
== bonds
[i
].aj()) ) ||
363 ( (ai
== bonds
[i
].aj()) && (aj
== bonds
[i
].ai()) ) )
365 bondlen
= bonds
[i
].c
; /* note: bonds[i].c might be NOTSET */
371 static real
get_angle(int nrang
, t_mybonded angles
[],
372 t_iatom ai
, t_iatom aj
, t_iatom ak
)
378 for (i
= 0; i
< nrang
&& (angle
== NOTSET
); i
++)
380 /* check both ways */
381 if ( ( (ai
== angles
[i
].ai()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ak()) ) ||
382 ( (ai
== angles
[i
].ak()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ai()) ) )
384 angle
= DEG2RAD
*angles
[i
].c
;
390 static char *get_atomtype_name_AB(t_atom
*atom
, gpp_atomtype_t atype
)
394 name
= get_atomtype_name(atom
->type
, atype
);
396 /* When using the decoupling option, atom types are changed
397 * to decoupled for the non-bonded interactions, but the virtual
398 * sites constructions should be based on the "normal" interactions.
399 * So we return the state B atom type name if the state A atom
400 * type is the decoupled one. We should actually check for the atom
401 * type number, but that's not passed here. So we check for
402 * the decoupled atom type name. This should not cause trouble
403 * as this code is only used for topologies with v-sites without
404 * parameters generated by pdb2gmx.
406 if (strcmp(name
, "decoupled") == 0)
408 name
= get_atomtype_name(atom
->typeB
, atype
);
414 static gmx_bool
calc_vsite3_param(gpp_atomtype_t atype
,
415 t_param
*param
, t_atoms
*at
,
416 int nrbond
, t_mybonded
*bonds
,
417 int nrang
, t_mybonded
*angles
)
419 /* i = virtual site | ,k
420 * j = 1st bonded heavy atom | i-j
421 * k,l = 2nd bonded atoms | `l
424 gmx_bool bXH3
, bError
;
425 real bjk
, bjl
, a
= -1, b
= -1;
426 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
427 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
431 for (i
= 0; i
< 4; i
++)
433 fprintf(debug
, "atom %d type %s ",
435 get_atomtype_name_AB(&at
->atom
[param
->a
[i
]], atype
));
437 fprintf(debug
, "\n");
440 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MNH", 3) == 0) &&
441 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MNH", 3) == 0) ) ||
442 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MCH3", 4) == 0) &&
443 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MCH3", 4) == 0) );
445 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
446 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
447 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
450 /* now we get some XH2/XH3 group specific construction */
451 /* note: we call the heavy atom 'C' and the X atom 'N' */
452 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, dM
, rM
;
455 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
456 bError
= bError
|| (bjk
!= bjl
);
458 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
459 aN
= std::max(param
->ak(), param
->al())+1;
461 /* get common bonds */
462 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
464 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
465 bError
= bError
|| (bMM
== NOTSET
) || (bCN
== NOTSET
);
467 /* calculate common things */
469 dM
= std::sqrt( sqr(bCM
) - sqr(rM
) );
471 /* are we dealing with the X atom? */
472 if (param
->ai() == aN
)
474 /* this is trivial */
475 a
= b
= 0.5 * bCN
/dM
;
480 /* get other bondlengths and angles: */
481 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
482 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
483 bError
= bError
|| (bNH
== NOTSET
) || (aCNH
== NOTSET
);
486 dH
= bCN
- bNH
* cos(aCNH
);
487 rH
= bNH
* sin(aCNH
);
489 a
= 0.5 * ( dH
/dM
+ rH
/rM
);
490 b
= 0.5 * ( dH
/dM
- rH
/rM
);
495 gmx_fatal(FARGS
, "calc_vsite3_param not implemented for the general case "
496 "(atom %d)", param
->ai()+1);
504 fprintf(debug
, "params for vsite3 %d: %g %g\n",
505 param
->ai()+1, param
->c0(), param
->c1());
511 static gmx_bool
calc_vsite3fd_param(t_param
*param
,
512 int nrbond
, t_mybonded
*bonds
,
513 int nrang
, t_mybonded
*angles
)
515 /* i = virtual site | ,k
516 * j = 1st bonded heavy atom | i-j
517 * k,l = 2nd bonded atoms | `l
521 real bij
, bjk
, bjl
, aijk
, aijl
, rk
, rl
;
523 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
524 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
525 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
526 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
527 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
528 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) ||
529 (aijk
== NOTSET
) || (aijl
== NOTSET
);
531 rk
= bjk
* sin(aijk
);
532 rl
= bjl
* sin(aijl
);
533 param
->c0() = rk
/ (rk
+ rl
);
534 param
->c1() = -bij
; /* 'bond'-length for fixed distance vsite */
538 fprintf(debug
, "params for vsite3fd %d: %g %g\n",
539 param
->ai()+1, param
->c0(), param
->c1());
544 static gmx_bool
calc_vsite3fad_param(t_param
*param
,
545 int nrbond
, t_mybonded
*bonds
,
546 int nrang
, t_mybonded
*angles
)
548 /* i = virtual site |
549 * j = 1st bonded heavy atom | i-j
550 * k = 2nd bonded heavy atom | `k-l
551 * l = 3d bonded heavy atom |
554 gmx_bool bSwapParity
, bError
;
557 bSwapParity
= ( param
->c1() == -1 );
559 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
560 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
561 bError
= (bij
== NOTSET
) || (aijk
== NOTSET
);
563 param
->c1() = bij
; /* 'bond'-length for fixed distance vsite */
564 param
->c0() = RAD2DEG
*aijk
; /* 'bond'-angle for fixed angle vsite */
568 param
->c0() = 360 - param
->c0();
573 fprintf(debug
, "params for vsite3fad %d: %g %g\n",
574 param
->ai()+1, param
->c0(), param
->c1());
579 static gmx_bool
calc_vsite3out_param(gpp_atomtype_t atype
,
580 t_param
*param
, t_atoms
*at
,
581 int nrbond
, t_mybonded
*bonds
,
582 int nrang
, t_mybonded
*angles
)
584 /* i = virtual site | ,k
585 * j = 1st bonded heavy atom | i-j
586 * k,l = 2nd bonded atoms | `l
587 * NOTE: i is out of the j-k-l plane!
590 gmx_bool bXH3
, bError
, bSwapParity
;
591 real bij
, bjk
, bjl
, aijk
, aijl
, akjl
, pijk
, pijl
, a
, b
, c
;
593 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
594 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
598 for (i
= 0; i
< 4; i
++)
600 fprintf(debug
, "atom %d type %s ",
601 param
->a
[i
]+1, get_atomtype_name_AB(&at
->atom
[param
->a
[i
]], atype
));
603 fprintf(debug
, "\n");
606 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MNH", 3) == 0) &&
607 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MNH", 3) == 0) ) ||
608 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MCH3", 4) == 0) &&
609 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MCH3", 4) == 0) );
611 /* check if construction parity must be swapped */
612 bSwapParity
= ( param
->c1() == -1 );
614 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
615 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
616 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
619 /* now we get some XH3 group specific construction */
620 /* note: we call the heavy atom 'C' and the X atom 'N' */
621 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, rHx
, rHy
, dM
, rM
;
624 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
625 bError
= bError
|| (bjk
!= bjl
);
627 /* the X atom (C or N) in the XH3 group is the first after the masses: */
628 aN
= std::max(param
->ak(), param
->al())+1;
630 /* get all bondlengths and angles: */
631 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
633 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
634 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
635 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
637 (bMM
== NOTSET
) || (bCN
== NOTSET
) || (bNH
== NOTSET
) || (aCNH
== NOTSET
);
640 dH
= bCN
- bNH
* cos(aCNH
);
641 rH
= bNH
* sin(aCNH
);
642 /* we assume the H's are symmetrically distributed */
643 rHx
= rH
*cos(DEG2RAD
*30);
644 rHy
= rH
*sin(DEG2RAD
*30);
646 dM
= std::sqrt( sqr(bCM
) - sqr(rM
) );
647 a
= 0.5*( (dH
/dM
) - (rHy
/rM
) );
648 b
= 0.5*( (dH
/dM
) + (rHy
/rM
) );
654 /* this is the general construction */
656 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
657 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
658 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
659 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
661 (bij
== NOTSET
) || (aijk
== NOTSET
) || (aijl
== NOTSET
) || (akjl
== NOTSET
);
663 pijk
= cos(aijk
)*bij
;
664 pijl
= cos(aijl
)*bij
;
665 a
= ( pijk
+ (pijk
*cos(akjl
)-pijl
) * cos(akjl
) / sqr(sin(akjl
)) ) / bjk
;
666 b
= ( pijl
+ (pijl
*cos(akjl
)-pijk
) * cos(akjl
) / sqr(sin(akjl
)) ) / bjl
;
667 c
= -std::sqrt( sqr(bij
) -
668 ( sqr(pijk
) - 2*pijk
*pijl
*cos(akjl
) + sqr(pijl
) )
670 / ( bjk
*bjl
*sin(akjl
) );
685 fprintf(debug
, "params for vsite3out %d: %g %g %g\n",
686 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
691 static gmx_bool
calc_vsite4fd_param(t_param
*param
,
692 int nrbond
, t_mybonded
*bonds
,
693 int nrang
, t_mybonded
*angles
)
695 /* i = virtual site | ,k
696 * j = 1st bonded heavy atom | i-j-m
697 * k,l,m = 2nd bonded atoms | `l
701 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
, akjm
, akjl
;
702 real pk
, pl
, pm
, cosakl
, cosakm
, sinakl
, sinakm
, cl
, cm
;
704 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
705 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
706 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
707 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
708 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
709 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
710 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
711 akjm
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->am());
712 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
713 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
714 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
) || (akjm
== NOTSET
) ||
722 cosakl
= (cos(akjl
) - cos(aijk
)*cos(aijl
)) / (sin(aijk
)*sin(aijl
));
723 cosakm
= (cos(akjm
) - cos(aijk
)*cos(aijm
)) / (sin(aijk
)*sin(aijm
));
724 if (cosakl
< -1 || cosakl
> 1 || cosakm
< -1 || cosakm
> 1)
726 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
727 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
728 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fd for atom %d: "
729 "cosakl=%f, cosakm=%f\n", param
->ai()+1, cosakl
, cosakm
);
731 sinakl
= std::sqrt(1-sqr(cosakl
));
732 sinakm
= std::sqrt(1-sqr(cosakm
));
734 /* note: there is a '+' because of the way the sines are calculated */
735 cl
= -pk
/ ( pl
*cosakl
- pk
+ pl
*sinakl
*(pm
*cosakm
-pk
)/(pm
*sinakm
) );
736 cm
= -pk
/ ( pm
*cosakm
- pk
+ pm
*sinakm
*(pl
*cosakl
-pk
)/(pl
*sinakl
) );
743 fprintf(debug
, "params for vsite4fd %d: %g %g %g\n",
744 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
753 calc_vsite4fdn_param(t_param
*param
,
754 int nrbond
, t_mybonded
*bonds
,
755 int nrang
, t_mybonded
*angles
)
757 /* i = virtual site | ,k
758 * j = 1st bonded heavy atom | i-j-m
759 * k,l,m = 2nd bonded atoms | `l
763 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
;
764 real pk
, pl
, pm
, a
, b
;
766 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
767 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
768 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
769 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
770 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
771 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
772 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
774 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
775 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
);
780 /* Calculate component of bond j-k along the direction i-j */
783 /* Calculate component of bond j-l along the direction i-j */
786 /* Calculate component of bond j-m along the direction i-j */
789 if (fabs(pl
) < 1000*GMX_REAL_MIN
|| fabs(pm
) < 1000*GMX_REAL_MIN
)
791 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
792 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
793 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fdn for atom %d: "
794 "pl=%f, pm=%f\n", param
->ai()+1, pl
, pm
);
806 fprintf(debug
, "params for vsite4fdn %d: %g %g %g\n",
807 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
816 int set_vsites(gmx_bool bVerbose
, t_atoms
*atoms
, gpp_atomtype_t atype
,
820 int nvsite
, nrbond
, nrang
, nridih
, nrset
;
821 gmx_bool bFirst
, bSet
, bERROR
;
822 at2vsitebond_t
*at2vb
;
831 fprintf(debug
, "\nCalculating parameters for virtual sites\n");
834 /* Make a reverse list to avoid ninteractions^2 operations */
835 at2vb
= make_at2vsitebond(atoms
->nr
, plist
);
837 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
839 if (interaction_function
[ftype
].flags
& IF_VSITE
)
841 nvsite
+= plist
[ftype
].nr
;
843 if (ftype
== F_VSITEN
)
845 /* We don't calculate parameters for VSITEN */
850 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
852 /* check if all parameters are set */
854 for (j
= 0; j
< NRFP(ftype
) && bSet
; j
++)
856 bSet
= plist
[ftype
].param
[i
].c
[j
] != NOTSET
;
861 fprintf(debug
, "bSet=%s ", bool_names
[bSet
]);
862 print_param(debug
, ftype
, i
, &plist
[ftype
].param
[i
]);
866 if (bVerbose
&& bFirst
)
868 fprintf(stderr
, "Calculating parameters for virtual sites\n");
872 nrbond
= nrang
= nridih
= 0;
877 /* now set the vsite parameters: */
878 get_bondeds(NRAL(ftype
), plist
[ftype
].param
[i
].a
, at2vb
,
879 &nrbond
, &bonds
, &nrang
, &angles
, &nridih
, &idihs
);
882 fprintf(debug
, "Found %d bonds, %d angles and %d idihs "
883 "for virtual site %d (%s)\n", nrbond
, nrang
, nridih
,
884 plist
[ftype
].param
[i
].ai()+1,
885 interaction_function
[ftype
].longname
);
886 print_bad(debug
, nrbond
, bonds
, nrang
, angles
, nridih
, idihs
);
892 calc_vsite3_param(atype
, &(plist
[ftype
].param
[i
]), atoms
,
893 nrbond
, bonds
, nrang
, angles
);
897 calc_vsite3fd_param(&(plist
[ftype
].param
[i
]),
898 nrbond
, bonds
, nrang
, angles
);
902 calc_vsite3fad_param(&(plist
[ftype
].param
[i
]),
903 nrbond
, bonds
, nrang
, angles
);
907 calc_vsite3out_param(atype
, &(plist
[ftype
].param
[i
]), atoms
,
908 nrbond
, bonds
, nrang
, angles
);
912 calc_vsite4fd_param(&(plist
[ftype
].param
[i
]),
913 nrbond
, bonds
, nrang
, angles
);
917 calc_vsite4fdn_param(&(plist
[ftype
].param
[i
]),
918 nrbond
, bonds
, nrang
, angles
);
921 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
923 interaction_function
[ftype
].longname
,
924 plist
[ftype
].param
[i
].ai()+1);
929 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
930 "for %s atom %d for this bonding configuration",
931 interaction_function
[ftype
].longname
,
932 plist
[ftype
].param
[i
].ai()+1);
939 if (debug
&& plist
[ftype
].nr
)
941 fprintf(stderr
, "Calculated parameters for %d out of %d %s atoms\n",
942 nrset
, plist
[ftype
].nr
, interaction_function
[ftype
].longname
);
947 done_at2vsitebond(atoms
->nr
, at2vb
);
952 void set_vsites_ptype(gmx_bool bVerbose
, gmx_moltype_t
*molt
)
961 fprintf(stderr
, "Setting particle type to V for virtual sites\n");
965 fprintf(stderr
, "checking %d functypes\n", F_NRE
);
967 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
969 il
= &molt
->ilist
[ftype
];
970 if (interaction_function
[ftype
].flags
& IF_VSITE
)
972 nra
= interaction_function
[ftype
].nratoms
;
978 fprintf(stderr
, "doing %d %s virtual sites\n",
979 (nrd
/ (nra
+1)), interaction_function
[ftype
].longname
);
982 for (i
= 0; (i
< nrd
); )
984 /* The virtual site */
986 molt
->atoms
.atom
[avsite
].ptype
= eptVSite
;
1000 static void check_vsite_constraints(t_params
*plist
,
1001 int cftype
, int vsite_type
[])
1008 ps
= &(plist
[cftype
]);
1009 for (i
= 0; (i
< ps
->nr
); i
++)
1011 for (k
= 0; k
< 2; k
++)
1013 atom
= ps
->param
[i
].a
[k
];
1014 if (vsite_type
[atom
] != NOTSET
)
1016 fprintf(stderr
, "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)\n",
1017 ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1, atom
+1);
1024 gmx_fatal(FARGS
, "There were %d virtual sites involved in constraints", n
);
1028 static void clean_vsite_bonds(t_params
*plist
, t_pindex pindex
[],
1029 int cftype
, int vsite_type
[])
1031 int ftype
, i
, j
, k
, m
, n
, nvsite
, nOut
, kept_i
;
1032 int nconverted
, nremoved
;
1033 atom_id atom
, oatom
, at1
, at2
;
1034 gmx_bool bKeep
, bRemove
, bUsed
, bPresent
, bThisFD
, bThisOUT
, bAllFD
, bFirstTwo
;
1037 if (cftype
== F_CONNBONDS
)
1042 ps
= &(plist
[cftype
]);
1047 for (i
= 0; (i
< ps
->nr
); i
++) /* for all bonds in the plist */
1050 const atom_id
*first_atoms
= NULL
;
1055 /* check if all virtual sites are constructed from the same atoms */
1059 fprintf(debug
, "constr %d %d:", ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1);
1061 for (k
= 0; (k
< 2) && !bKeep
&& !bRemove
; k
++)
1063 /* for all atoms in the bond */
1064 atom
= ps
->param
[i
].a
[k
];
1065 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1068 bThisFD
= ( (pindex
[atom
].ftype
== F_VSITE3FD
) ||
1069 (pindex
[atom
].ftype
== F_VSITE3FAD
) ||
1070 (pindex
[atom
].ftype
== F_VSITE4FD
) ||
1071 (pindex
[atom
].ftype
== F_VSITE4FDN
) );
1072 bThisOUT
= ( (pindex
[atom
].ftype
== F_VSITE3OUT
) &&
1073 (interaction_function
[cftype
].flags
& IF_CONSTRAINT
) );
1074 bAllFD
= bAllFD
&& bThisFD
;
1075 if (bThisFD
|| bThisOUT
)
1079 fprintf(debug
, " %s", bThisOUT
? "out" : "fd");
1081 oatom
= ps
->param
[i
].a
[1-k
]; /* the other atom */
1082 if (vsite_type
[oatom
] == NOTSET
&&
1083 vsite_type
[oatom
] != F_VSITEN
&&
1084 oatom
== plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].aj())
1086 /* if the other atom isn't a vsite, and it is AI */
1094 fprintf(debug
, " D-AI");
1100 /* TODO This fragment, and corresponding logic in
1101 clean_vsite_angles and clean_vsite_dihs should
1102 be refactored into a common function */
1105 /* if this is the first vsite we encounter then
1106 store construction atoms */
1107 /* TODO This would be nicer to implement with
1108 a C++ "vector view" class" with an
1109 STL-container-like interface. */
1110 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1111 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1115 GMX_ASSERT(vsnral
!= 0, "nvsite > 1 must have vsnral != 0");
1116 GMX_ASSERT(first_atoms
!= NULL
, "nvsite > 1 must have first_atoms != NULL");
1117 /* if it is not the first then
1118 check if this vsite is constructed from the same atoms */
1119 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1121 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1123 const atom_id
*atoms
;
1126 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1127 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1129 if (atoms
[m
] == first_atoms
[n
])
1139 fprintf(debug
, " !present");
1149 fprintf(debug
, " !same#at");
1163 /* if we have no virtual sites in this bond, keep it */
1168 fprintf(debug
, " no vsite");
1173 /* TODO This loop and the corresponding loop in
1174 check_vsite_angles should be refactored into a common
1176 /* check if all non-vsite atoms are used in construction: */
1178 for (k
= 0; (k
< 2) && !bKeep
; k
++) /* for all atoms in the bond */
1180 atom
= ps
->param
[i
].a
[k
];
1181 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1184 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1186 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1188 if (atom
== first_atoms
[m
])
1191 bFirstTwo
= bFirstTwo
&& m
< 2;
1199 fprintf(debug
, " !used");
1205 if (!( bAllFD
&& bFirstTwo
) )
1207 /* Two atom bonded interactions include constraints.
1208 * We need to remove constraints between vsite pairs that have
1209 * a fixed distance due to being constructed from the same
1210 * atoms, since this can be numerically unstable.
1212 for (m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1214 at1
= first_atoms
[m
];
1215 at2
= first_atoms
[(m
+1) % vsnral
];
1217 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1219 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1221 for (j
= 0; (j
< plist
[ftype
].nr
) && !bPresent
; j
++)
1223 /* all constraints until one matches */
1224 bPresent
= ( ( (plist
[ftype
].param
[j
].ai() == at1
) &&
1225 (plist
[ftype
].param
[j
].aj() == at2
) ) ||
1226 ( (plist
[ftype
].param
[j
].ai() == at2
) &&
1227 (plist
[ftype
].param
[j
].aj() == at1
) ) );
1236 fprintf(debug
, " !bonded");
1247 fprintf(debug
, " keeping");
1249 /* now copy the bond to the new array */
1250 ps
->param
[kept_i
] = ps
->param
[i
];
1253 else if (IS_CHEMBOND(cftype
))
1255 srenew(plist
[F_CONNBONDS
].param
, plist
[F_CONNBONDS
].nr
+1);
1256 plist
[F_CONNBONDS
].param
[plist
[F_CONNBONDS
].nr
] = ps
->param
[i
];
1257 plist
[F_CONNBONDS
].nr
++;
1266 fprintf(debug
, "\n");
1272 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1273 nremoved
, interaction_function
[cftype
].longname
, kept_i
);
1277 fprintf(stderr
, "Converted %4d %15ss with virtual sites to connections, %5d left\n",
1278 nconverted
, interaction_function
[cftype
].longname
, kept_i
);
1282 fprintf(stderr
, "Warning: removed %d %ss with vsite with %s construction\n"
1283 " This vsite construction does not guarantee constant "
1285 " If the constructions were generated by pdb2gmx ignore "
1287 nOut
, interaction_function
[cftype
].longname
,
1288 interaction_function
[F_VSITE3OUT
].longname
);
1293 static void clean_vsite_angles(t_params
*plist
, t_pindex pindex
[],
1294 int cftype
, int vsite_type
[],
1295 at2vsitecon_t
*at2vc
)
1297 int i
, j
, k
, m
, n
, nvsite
, kept_i
;
1298 atom_id atom
, at1
, at2
;
1299 gmx_bool bKeep
, bUsed
, bPresent
, bAll3FAD
, bFirstTwo
;
1302 ps
= &(plist
[cftype
]);
1304 for (i
= 0; (i
< ps
->nr
); i
++) /* for all angles in the plist */
1307 const atom_id
*first_atoms
= NULL
;
1311 /* check if all virtual sites are constructed from the same atoms */
1313 for (k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1315 atom
= ps
->param
[i
].a
[k
];
1316 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1319 bAll3FAD
= bAll3FAD
&& (pindex
[atom
].ftype
== F_VSITE3FAD
);
1322 /* store construction atoms of first vsite */
1323 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1324 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1328 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1329 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1330 /* check if this vsite is constructed from the same atoms */
1331 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1333 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1335 const atom_id
*atoms
;
1338 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1339 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1341 if (atoms
[m
] == first_atoms
[n
])
1360 /* keep all angles with no virtual sites in them or
1361 with virtual sites with more than 3 constr. atoms */
1362 if (nvsite
== 0 && vsnral
> 3)
1367 /* check if all non-vsite atoms are used in construction: */
1369 for (k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1371 atom
= ps
->param
[i
].a
[k
];
1372 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1375 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1377 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1379 if (atom
== first_atoms
[m
])
1382 bFirstTwo
= bFirstTwo
&& m
< 2;
1392 if (!( bAll3FAD
&& bFirstTwo
) )
1394 /* check if all constructing atoms are constrained together */
1395 for (m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1397 at1
= first_atoms
[m
];
1398 at2
= first_atoms
[(m
+1) % vsnral
];
1400 for (j
= 0; j
< at2vc
[at1
].nr
; j
++)
1402 if (at2vc
[at1
].aj
[j
] == at2
)
1416 /* now copy the angle to the new array */
1417 ps
->param
[kept_i
] = ps
->param
[i
];
1422 if (ps
->nr
!= kept_i
)
1424 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1425 ps
->nr
-kept_i
, interaction_function
[cftype
].longname
, kept_i
);
1430 static void clean_vsite_dihs(t_params
*plist
, t_pindex pindex
[],
1431 int cftype
, int vsite_type
[])
1436 ps
= &(plist
[cftype
]);
1439 for (i
= 0; (i
< ps
->nr
); i
++) /* for all dihedrals in the plist */
1441 int k
, m
, n
, nvsite
;
1443 const atom_id
*first_atoms
= NULL
;
1445 gmx_bool bKeep
, bUsed
, bPresent
;
1449 /* check if all virtual sites are constructed from the same atoms */
1451 for (k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1453 atom
= ps
->param
[i
].a
[k
];
1454 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1458 /* store construction atoms of first vsite */
1459 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1460 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1463 fprintf(debug
, "dih w. vsite: %d %d %d %d\n",
1464 ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1,
1465 ps
->param
[i
].ak()+1, ps
->param
[i
].al()+1);
1466 fprintf(debug
, "vsite %d from: %d %d %d\n",
1467 atom
+1, first_atoms
[0]+1, first_atoms
[1]+1, first_atoms
[2]+1);
1472 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1473 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1474 /* check if this vsite is constructed from the same atoms */
1475 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1477 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1479 const atom_id
*atoms
;
1482 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1483 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1485 if (atoms
[m
] == first_atoms
[n
])
1497 /* TODO clean_site_bonds and _angles do this increment
1498 at the top of the loop. Refactor this for
1504 /* keep all dihedrals with no virtual sites in them */
1510 /* check if all atoms in dihedral are either virtual sites, or used in
1511 construction of virtual sites. If so, keep it, if not throw away: */
1512 for (k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1514 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1515 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1516 atom
= ps
->param
[i
].a
[k
];
1517 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1519 /* vsnral will be set here, we don't get here with nvsite==0 */
1521 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1523 if (atom
== first_atoms
[m
])
1533 fprintf(debug
, "unused atom in dih: %d\n", atom
+1);
1541 ps
->param
[kept_i
] = ps
->param
[i
];
1546 if (ps
->nr
!= kept_i
)
1548 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1549 ps
->nr
-kept_i
, interaction_function
[cftype
].longname
, kept_i
);
1554 void clean_vsite_bondeds(t_params
*plist
, int natoms
, gmx_bool bRmVSiteBds
)
1556 int i
, k
, nvsite
, ftype
, vsite
, parnr
;
1559 at2vsitecon_t
*at2vc
;
1561 pindex
= 0; /* avoid warnings */
1562 /* make vsite_type array */
1563 snew(vsite_type
, natoms
);
1564 for (i
= 0; i
< natoms
; i
++)
1566 vsite_type
[i
] = NOTSET
;
1569 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1571 if (interaction_function
[ftype
].flags
& IF_VSITE
)
1573 nvsite
+= plist
[ftype
].nr
;
1575 while (i
< plist
[ftype
].nr
)
1577 vsite
= plist
[ftype
].param
[i
].ai();
1578 if (vsite_type
[vsite
] == NOTSET
)
1580 vsite_type
[vsite
] = ftype
;
1584 gmx_fatal(FARGS
, "multiple vsite constructions for atom %d", vsite
+1);
1586 if (ftype
== F_VSITEN
)
1588 while (i
< plist
[ftype
].nr
&& plist
[ftype
].param
[i
].ai() == vsite
)
1601 /* the rest only if we have virtual sites: */
1604 fprintf(stderr
, "Cleaning up constraints %swith virtual sites\n",
1605 bRmVSiteBds
? "and constant bonded interactions " : "");
1607 /* Make a reverse list to avoid ninteractions^2 operations */
1608 at2vc
= make_at2vsitecon(natoms
, plist
);
1610 snew(pindex
, natoms
);
1611 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1613 /* Here we skip VSITEN. In neary all practical use cases this
1614 * is not an issue, since VSITEN is intended for constructing
1615 * additional interaction sites, not for replacing normal atoms
1616 * with bonded interactions. Thus we do not expect constant
1617 * bonded interactions. If VSITEN does get used with constant
1618 * bonded interactions, these are not removed which only leads
1619 * to very minor extra computation and constant energy.
1620 * The only problematic case is onstraints between VSITEN
1621 * constructions with fixed distance (which is anyhow useless).
1622 * This will generate a fatal error in check_vsite_constraints.
1624 if ((interaction_function
[ftype
].flags
& IF_VSITE
) &&
1627 for (parnr
= 0; (parnr
< plist
[ftype
].nr
); parnr
++)
1629 k
= plist
[ftype
].param
[parnr
].ai();
1630 pindex
[k
].ftype
= ftype
;
1631 pindex
[k
].parnr
= parnr
;
1638 for (i
= 0; i
< natoms
; i
++)
1640 fprintf(debug
, "atom %d vsite_type %s\n", i
,
1641 vsite_type
[i
] == NOTSET
? "NOTSET" :
1642 interaction_function
[vsite_type
[i
]].name
);
1646 /* remove interactions that include virtual sites */
1647 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1649 if ( ( ( interaction_function
[ftype
].flags
& IF_BOND
) && bRmVSiteBds
) ||
1650 ( interaction_function
[ftype
].flags
& IF_CONSTRAINT
) )
1652 if (interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
) )
1654 clean_vsite_bonds (plist
, pindex
, ftype
, vsite_type
);
1656 else if (interaction_function
[ftype
].flags
& IF_ATYPE
)
1658 clean_vsite_angles(plist
, pindex
, ftype
, vsite_type
, at2vc
);
1660 else if ( (ftype
== F_PDIHS
) || (ftype
== F_IDIHS
) )
1662 clean_vsite_dihs (plist
, pindex
, ftype
, vsite_type
);
1666 /* check that no remaining constraints include virtual sites */
1667 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1669 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1671 check_vsite_constraints(plist
, ftype
, vsite_type
);
1675 done_at2vsitecon(natoms
, at2vc
);