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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /* This file is completely threadsafe - keep it that way! */
55 #include "gmx_fatal.h"
57 static void copy_bond (t_params *pr, int to, int from)
58 /* copies an entry in a bond list to another position.
59 * does no allocing or freeing of memory
62 /*memcpy((char*) &(pr->param[to]),(char*) &(pr->param[from]),
63 (size_t)sizeof(pr->param[from]));*/
67 range_check(to,0,pr->nr);
68 range_check(from,0,pr->nr);
69 for(i=0; (i<MAXATOMLIST); i++)
70 pr->param[to].a[i] = pr->param[from].a[i];
71 for(i=0; (i<MAXFORCEPARAM); i++)
72 pr->param[to].c[i] = pr->param[from].c[i];
73 for(i=0; (i<MAXSLEN); i++)
74 pr->param[to].s[i] = pr->param[from].s[i];
78 static int count_hydrogens (char ***atomname, int nra, atom_id a[])
83 gmx_fatal(FARGS,"Cannot call count_hydrogens with no atomname (%s %d)",
87 for (i=0; (i<nra); i++)
88 if (toupper(**(atomname[a[i]]))=='H')
94 void make_shake (t_params plist[],t_atoms *atoms,gpp_atomtype_t at,int nshake)
96 char ***info=atoms->atomname;
101 int nb,b,i,j,ftype,ftype_a;
104 if (nshake != eshNONE) {
107 printf("turning H bonds into constraints...\n");
110 printf("turning all bonds into constraints...\n");
113 printf("turning all bonds and H angles into constraints...\n");
116 printf("turning all bonds and angles into constraints...\n");
119 gmx_fatal(FARGS,"Invalid option for make_shake (%d)",nshake);
122 if ((nshake == eshHANGLES) || (nshake == eshALLANGLES)) {
123 /* Add all the angles with hydrogens to the shake list
124 * and remove them from the bond list
126 for(ftype = 0; (ftype < F_NRE); ftype++) {
127 if (interaction_function[ftype].flags & IF_BTYPE) {
128 bonds=&(plist[ftype]);
130 for(ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) {
131 if (interaction_function[ftype_a].flags & IF_ATYPE) {
132 pr = &(plist[ftype_a]);
134 for (i=0; (i < pr->nr); ) {
139 printf("Angle: %d-%d-%d\n",ang->AI,ang->AJ,ang->AK);
141 numhydrogens = count_hydrogens(info,3,ang->a);
142 if ((nshake == eshALLANGLES) ||
143 (numhydrogens > 1) ||
144 (numhydrogens == 1 && toupper(**(info[ang->a[1]]))=='O')) {
145 /* Can only add hydrogen angle shake, if the two bonds
147 * append this angle to the shake list
152 /* Calculate length of constraint */
155 for (j=0; !bFound && (j<bonds->nr); j++) {
156 bond=&(bonds->param[j]);
157 if (((bond->AI==ang->AI) &&
158 (bond->AJ==ang->AJ)) ||
159 ((bond->AI==ang->AJ) &&
160 (bond->AJ==ang->AI)))
162 if (((bond->AI==ang->AK) &&
163 (bond->AJ==ang->AJ)) ||
164 ((bond->AI==ang->AJ) &&
165 (bond->AJ==ang->AK)))
167 bFound = (b_ij!=0.0) && (b_jk!=0.0);
170 /* apply law of cosines */
171 p.C0 = sqrt( b_ij*b_ij + b_jk*b_jk -
172 2.0*b_ij*b_jk*cos(DEG2RAD*ang->C0) );
175 printf("p: %d, q: %d, dist: %12.5e\n",p.AI,p.AJ,p.C0);
177 add_param_to_list (&(plist[F_CONSTR]),&p);
178 /* move the last bond to this position */
179 copy_bond (pr,i,pr->nr-1);
180 /* should free memory here!! */
191 } /* if shake angles */
193 /* Add all the bonds with hydrogens to the shake list
194 * and remove them from the bond list
196 for (ftype=0; (ftype < F_NRE); ftype++) {
197 if (interaction_function[ftype].flags & IF_BTYPE) {
198 pr = &(plist[ftype]);
200 for(i=0; i<pr->nr; i++) {
201 if ( (nshake != eshHBONDS) ||
202 (count_hydrogens (info,2,pr->param[i].a) > 0) ) {
203 /* append this bond to the shake list */
204 p.AI = pr->param[i].AI;
205 p.AJ = pr->param[i].AJ;
206 p.C0 = pr->param[i].C0;
207 p.C1 = pr->param[i].C2;
208 add_param_to_list (&(plist[F_CONSTR]),&p);