3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35 /* This file is completely threadsafe - keep it that way! */
52 #include "gmx_fatal.h"
54 static void copy_bond (t_params *pr, int to, int from)
55 /* copies an entry in a bond list to another position.
56 * does no allocing or freeing of memory
59 /*memcpy((char*) &(pr->param[to]),(char*) &(pr->param[from]),
60 (size_t)sizeof(pr->param[from]));*/
64 range_check(to,0,pr->nr);
65 range_check(from,0,pr->nr);
66 for(i=0; (i<MAXATOMLIST); i++)
67 pr->param[to].a[i] = pr->param[from].a[i];
68 for(i=0; (i<MAXFORCEPARAM); i++)
69 pr->param[to].c[i] = pr->param[from].c[i];
70 for(i=0; (i<MAXSLEN); i++)
71 pr->param[to].s[i] = pr->param[from].s[i];
75 static int count_hydrogens (char ***atomname, int nra, atom_id a[])
80 gmx_fatal(FARGS,"Cannot call count_hydrogens with no atomname (%s %d)",
84 for (i=0; (i<nra); i++)
85 if (toupper(**(atomname[a[i]]))=='H')
91 void make_shake (t_params plist[],t_atoms *atoms,gpp_atomtype_t at,int nshake)
93 char ***info=atoms->atomname;
98 int nb,b,i,j,ftype,ftype_a;
101 if (nshake != eshNONE) {
104 printf("turning H bonds into constraints...\n");
107 printf("turning all bonds into constraints...\n");
110 printf("turning all bonds and H angles into constraints...\n");
113 printf("turning all bonds and angles into constraints...\n");
116 gmx_fatal(FARGS,"Invalid option for make_shake (%d)",nshake);
119 if ((nshake == eshHANGLES) || (nshake == eshALLANGLES)) {
120 /* Add all the angles with hydrogens to the shake list
121 * and remove them from the bond list
123 for(ftype = 0; (ftype < F_NRE); ftype++) {
124 if (interaction_function[ftype].flags & IF_BTYPE) {
125 bonds=&(plist[ftype]);
127 for(ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) {
128 if (interaction_function[ftype_a].flags & IF_ATYPE) {
129 pr = &(plist[ftype_a]);
131 for (i=0; (i < pr->nr); ) {
136 printf("Angle: %d-%d-%d\n",ang->AI,ang->AJ,ang->AK);
138 numhydrogens = count_hydrogens(info,3,ang->a);
139 if ((nshake == eshALLANGLES) ||
140 (numhydrogens > 1) ||
141 (numhydrogens == 1 && toupper(**(info[ang->a[1]]))=='O')) {
142 /* Can only add hydrogen angle shake, if the two bonds
144 * append this angle to the shake list
149 /* Calculate length of constraint */
152 for (j=0; !bFound && (j<bonds->nr); j++) {
153 bond=&(bonds->param[j]);
154 if (((bond->AI==ang->AI) &&
155 (bond->AJ==ang->AJ)) ||
156 ((bond->AI==ang->AJ) &&
157 (bond->AJ==ang->AI)))
159 if (((bond->AI==ang->AK) &&
160 (bond->AJ==ang->AJ)) ||
161 ((bond->AI==ang->AJ) &&
162 (bond->AJ==ang->AK)))
164 bFound = (b_ij!=0.0) && (b_jk!=0.0);
167 /* apply law of cosines */
168 p.C0 = sqrt( b_ij*b_ij + b_jk*b_jk -
169 2.0*b_ij*b_jk*cos(DEG2RAD*ang->C0) );
172 printf("p: %d, q: %d, dist: %12.5e\n",p.AI,p.AJ,p.C0);
174 add_param_to_list (&(plist[F_CONSTR]),&p);
175 /* move the last bond to this position */
176 copy_bond (pr,i,pr->nr-1);
177 /* should free memory here!! */
188 } /* if shake angles */
190 /* Add all the bonds with hydrogens to the shake list
191 * and remove them from the bond list
193 for (ftype=0; (ftype < F_NRE); ftype++) {
194 if (interaction_function[ftype].flags & IF_BTYPE) {
195 pr = &(plist[ftype]);
197 for(i=0; i<pr->nr; i++) {
198 if ( (nshake != eshHBONDS) ||
199 (count_hydrogens (info,2,pr->param[i].a) > 0) ) {
200 /* append this bond to the shake list */
201 p.AI = pr->param[i].AI;
202 p.AJ = pr->param[i].AJ;
203 p.C0 = pr->param[i].C0;
204 p.C1 = pr->param[i].C2;
205 add_param_to_list (&(plist[F_CONSTR]),&p);