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,2013, 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.
55 #include "gmx_fatal.h"
57 #include "gpp_atomtype.h"
58 #include "gpp_bond_atomtype.h"
60 void generate_nbparams(int comb,int ftype,t_params *plist,gpp_atomtype_t atype,
65 real c,bi,bj,ci,cj,ci0,ci1,ci2,cj0,cj1,cj2;
68 /* Lean mean shortcuts */
69 nr = get_atomtype_ntypes(atype);
71 snew(plist->param,nr*nr);
74 /* Fill the matrix with force parameters */
80 for(i=k=0; (i<nr); i++) {
81 for(j=0; (j<nr); j++,k++) {
82 for(nf=0; (nf<nrfp); nf++) {
83 ci = get_atomtype_nbparam(i,nf,atype);
84 cj = get_atomtype_nbparam(j,nf,atype);
86 plist->param[k].c[nf] = c;
92 case eCOMB_ARITHMETIC:
93 /* c0 and c1 are epsilon and sigma */
94 for (i=k=0; (i < nr); i++)
95 for (j=0; (j < nr); j++,k++) {
96 ci0 = get_atomtype_nbparam(i,0,atype);
97 cj0 = get_atomtype_nbparam(j,0,atype);
98 ci1 = get_atomtype_nbparam(i,1,atype);
99 cj1 = get_atomtype_nbparam(j,1,atype);
100 plist->param[k].c[0] = (ci0+cj0)*0.5;
101 plist->param[k].c[1] = sqrt(ci1*cj1);
105 case eCOMB_GEOM_SIG_EPS:
106 /* c0 and c1 are epsilon and sigma */
107 for (i=k=0; (i < nr); i++)
108 for (j=0; (j < nr); j++,k++) {
109 ci0 = get_atomtype_nbparam(i,0,atype);
110 cj0 = get_atomtype_nbparam(j,0,atype);
111 ci1 = get_atomtype_nbparam(i,1,atype);
112 cj1 = get_atomtype_nbparam(j,1,atype);
113 plist->param[k].c[0] = sqrt(ci0*cj0);
114 plist->param[k].c[1] = sqrt(ci1*cj1);
119 gmx_fatal(FARGS,"No such combination rule %d",comb);
122 gmx_incons("Topology processing, generate nb parameters");
126 /* Buckingham rules */
127 for(i=k=0; (i<nr); i++)
128 for(j=0; (j<nr); j++,k++) {
129 ci0 = get_atomtype_nbparam(i,0,atype);
130 cj0 = get_atomtype_nbparam(j,0,atype);
131 ci2 = get_atomtype_nbparam(i,2,atype);
132 cj2 = get_atomtype_nbparam(j,2,atype);
133 bi = get_atomtype_nbparam(i,1,atype);
134 bj = get_atomtype_nbparam(j,1,atype);
135 plist->param[k].c[0] = sqrt(ci0 * cj0);
136 if ((bi == 0) || (bj == 0))
137 plist->param[k].c[1] = 0;
139 plist->param[k].c[1] = 2.0/(1/bi+1/bj);
140 plist->param[k].c[2] = sqrt(ci2 * cj2);
145 sprintf(errbuf,"Invalid nonbonded type %s",
146 interaction_function[ftype].longname);
147 warning_error(wi,errbuf);
151 static void realloc_nb_params(gpp_atomtype_t at,
152 t_nbparam ***nbparam,t_nbparam ***pair)
154 /* Add space in the non-bonded parameters matrix */
155 int atnr = get_atomtype_ntypes(at);
156 srenew(*nbparam,atnr);
157 snew((*nbparam)[atnr-1],atnr);
160 snew((*pair)[atnr-1],atnr);
164 static void copy_B_from_A(int ftype,double *c)
168 nrfpA = NRFPA(ftype);
169 nrfpB = NRFPB(ftype);
171 /* Copy the B parameters from the first nrfpB A parameters */
172 for(i=0; (i<nrfpB); i++) {
177 void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
178 char *line,int nb_funct,
179 t_nbparam ***nbparam,t_nbparam ***pair,
186 t_xlate xl[eptNR] = {
194 int nr,i,nfields,j,pt,nfp0=-1;
196 char type[STRLEN],btype[STRLEN],ptype[STRLEN];
198 double c[MAXFORCEPARAM];
199 double radius,vol,surftens,gb_radius,S_hct;
200 char tmpfield[12][100]; /* Max 12 fields of width 100 */
205 gmx_bool have_atomic_number;
206 gmx_bool have_bonded_type;
211 /* First assign input line to temporary array */
212 nfields=sscanf(line,"%s%s%s%s%s%s%s%s%s%s%s%s",
213 tmpfield[0],tmpfield[1],tmpfield[2],tmpfield[3],tmpfield[4],tmpfield[5],
214 tmpfield[6],tmpfield[7],tmpfield[8],tmpfield[9],tmpfield[10],tmpfield[11]);
216 /* Comments on optional fields in the atomtypes section:
218 * The force field format is getting a bit old. For OPLS-AA we needed
219 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
220 * we also needed the atomic numbers.
221 * To avoid making all old or user-generated force fields unusable we
222 * have introduced both these quantities as optional columns, and do some
223 * acrobatics to check whether they are present or not.
224 * This will all look much nicer when we switch to XML... sigh.
226 * Field 0 (mandatory) is the nonbonded type name. (string)
227 * Field 1 (optional) is the bonded type (string)
228 * Field 2 (optional) is the atomic number (int)
229 * Field 3 (mandatory) is the mass (numerical)
230 * Field 4 (mandatory) is the charge (numerical)
231 * Field 5 (mandatory) is the particle type (single character)
232 * This is followed by a number of nonbonded parameters.
234 * The safest way to identify the format is the particle type field.
236 * So, here is what we do:
238 * A. Read in the first six fields as strings
239 * B. If field 3 (starting from 0) is a single char, we have neither
240 * bonded_type or atomic numbers.
241 * C. If field 5 is a single char we have both.
242 * D. If field 4 is a single char we check field 1. If this begins with
243 * an alphabetical character we have bonded types, otherwise atomic numbers.
252 if( (strlen(tmpfield[5])==1) && isalpha(tmpfield[5][0]) )
254 have_bonded_type = TRUE;
255 have_atomic_number = TRUE;
257 else if( (strlen(tmpfield[3])==1) && isalpha(tmpfield[3][0]) )
259 have_bonded_type = FALSE;
260 have_atomic_number = FALSE;
264 have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 );
265 have_atomic_number = !have_bonded_type;
268 /* optional fields */
281 if( have_atomic_number )
283 if ( have_bonded_type )
285 nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
286 type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],
287 &radius,&vol,&surftens,&gb_radius);
296 /* have_atomic_number && !have_bonded_type */
297 nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
298 type,&atomnr,&m,&q,ptype,&c[0],&c[1],
299 &radius,&vol,&surftens,&gb_radius);
309 if ( have_bonded_type )
311 /* !have_atomic_number && have_bonded_type */
312 nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
313 type,btype,&m,&q,ptype,&c[0],&c[1],
314 &radius,&vol,&surftens,&gb_radius);
323 /* !have_atomic_number && !have_bonded_type */
324 nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
325 type,&m,&q,ptype,&c[0],&c[1],
326 &radius,&vol,&surftens,&gb_radius);
335 if( !have_bonded_type )
340 if( !have_atomic_number )
350 if( have_atomic_number )
352 if ( have_bonded_type )
354 nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
355 type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
356 &radius,&vol,&surftens,&gb_radius);
365 /* have_atomic_number && !have_bonded_type */
366 nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
367 type,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
368 &radius,&vol,&surftens,&gb_radius);
378 if ( have_bonded_type )
380 /* !have_atomic_number && have_bonded_type */
381 nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
382 type,btype,&m,&q,ptype,&c[0],&c[1],&c[2],
383 &radius,&vol,&surftens,&gb_radius);
392 /* !have_atomic_number && !have_bonded_type */
393 nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
394 type,&m,&q,ptype,&c[0],&c[1],&c[2],
395 &radius,&vol,&surftens,&gb_radius);
404 if( !have_bonded_type )
409 if( !have_atomic_number )
417 gmx_fatal(FARGS,"Invalid function type %d in push_at %s %d",nb_funct,
420 for(j=nfp0; (j<MAXFORCEPARAM); j++)
423 if(strlen(type)==1 && isdigit(type[0]))
424 gmx_fatal(FARGS,"Atom type names can't be single digits.");
426 if(strlen(btype)==1 && isdigit(btype[0]))
427 gmx_fatal(FARGS,"Bond atom type names can't be single digits.");
429 /* Hack to read old topologies */
430 if (gmx_strcasecmp(ptype,"D") == 0)
432 for(j=0; (j<eptNR); j++)
433 if (gmx_strcasecmp(ptype,xl[j].entry) == 0)
436 gmx_fatal(FARGS,"Invalid particle type %s on line %s",
440 fprintf(debug,"ptype: %s\n",ptype_str[pt]);
445 for (i=0; (i<MAXFORCEPARAM); i++)
448 if ((batype_nr = get_bond_atomtype_type(btype,bat)) == NOTSET)
449 add_bond_atomtype(bat,symtab,btype);
450 batype_nr = get_bond_atomtype_type(btype,bat);
452 if ((nr = get_atomtype_type(type,at)) != NOTSET) {
453 sprintf(errbuf,"Overriding atomtype %s",type);
455 if ((nr = set_atomtype(nr,at,symtab,atom,type,param,batype_nr,
456 radius,vol,surftens,atomnr,gb_radius,S_hct)) == NOTSET)
457 gmx_fatal(FARGS,"Replacing atomtype %s failed",type);
459 else if ((nr = add_atomtype(at,symtab,atom,type,param,
460 batype_nr,radius,vol,
461 surftens,atomnr,gb_radius,S_hct)) == NOTSET)
462 gmx_fatal(FARGS,"Adding atomtype %s failed",type);
464 /* Add space in the non-bonded parameters matrix */
465 realloc_nb_params(at,nbparam,pair);
471 static void push_bondtype(t_params * bt,
475 gmx_bool bAllowRepeat,
480 gmx_bool bTest,bFound,bCont,bId;
482 int nrfp = NRFP(ftype);
485 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
486 are on directly _adjacent_ lines.
489 /* First check if our atomtypes are _identical_ (not reversed) to the previous
490 entry. If they are not identical we search for earlier duplicates. If they are
491 we can skip it, since we already searched for the first line
498 if(bAllowRepeat && nr>1)
500 for (j=0,bCont=TRUE; (j < nral); j++)
502 bCont=bCont && (b->a[j] == bt->param[nr-2].a[j]);
506 /* Search for earlier duplicates if this entry was not a continuation
507 from the previous line.
512 for (i=0; (i < nr); i++) {
514 for (j=0; (j < nral); j++)
515 bTest=(bTest && (b->a[j] == bt->param[i].a[j]));
518 for(j=0; (j<nral); j++)
519 bTest=(bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
524 for (j=0; (j < nrfp); j++)
525 bId = bId && (bt->param[i].c[j] == b->c[j]);
527 sprintf(errbuf,"Overriding %s parameters.%s",
528 interaction_function[ftype].longname,
529 (ftype==F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
531 fprintf(stderr," old:");
532 for (j=0; (j < nrfp); j++)
533 fprintf(stderr," %g",bt->param[i].c[j]);
534 fprintf(stderr," \n new: %s\n\n",line);
538 for (j=0; (j < nrfp); j++)
539 bt->param[i].c[j] = b->c[j];
548 /* fill the arrays up and down */
549 memcpy(bt->param[bt->nr].c, b->c,sizeof(b->c));
550 memcpy(bt->param[bt->nr].a, b->a,sizeof(b->a));
551 memcpy(bt->param[bt->nr+1].c,b->c,sizeof(b->c));
553 /* The definitions of linear angles depend on the order of atoms,
554 * that means that for atoms i-j-k, with certain parameter a, the
555 * corresponding k-j-i angle will have parameter 1-a.
557 if (ftype == F_LINEAR_ANGLES)
559 bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
560 bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
563 for (j=0; (j < nral); j++)
565 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
572 void push_bt(directive d,t_params bt[],int nral,
574 t_bond_atomtype bat,char *line,
577 const char *formal[MAXATOMLIST+1] = {
586 const char *formnl[MAXATOMLIST+1] = {
592 "%*s%*s%*s%*s%*s%*s",
593 "%*s%*s%*s%*s%*s%*s%*s"
595 const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
596 int i,ft,ftype,nn,nrfp,nrfpA,nrfpB;
598 char alc[MAXATOMLIST+1][20];
599 /* One force parameter more, so we can check if we read too many */
600 double c[MAXFORCEPARAM+1];
604 if ((bat && at) || (!bat && !at))
605 gmx_incons("You should pass either bat or at to push_bt");
607 /* Make format string (nral ints+functype) */
608 if ((nn=sscanf(line,formal[nral],
609 alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
610 sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
611 warning_error(wi,errbuf);
615 ft = strtol(alc[nral],NULL,10);
616 ftype = ifunc_index(d,ft);
618 nrfpA = interaction_function[ftype].nrfpA;
619 nrfpB = interaction_function[ftype].nrfpB;
620 strcpy(f1,formnl[nral]);
622 if ((nn=sscanf(line,f1,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5],&c[6],&c[7],&c[8],&c[9],&c[10],&c[11],&c[12]))
625 /* Copy the B-state from the A-state */
626 copy_B_from_A(ftype,c);
629 warning_error(wi,"Not enough parameters");
630 } else if (nn > nrfpA && nn < nrfp) {
631 warning_error(wi,"Too many parameters or not enough parameters for topology B");
632 } else if (nn > nrfp) {
633 warning_error(wi,"Too many parameters");
635 for(i=nn; (i<nrfp); i++)
639 for(i=0; (i<nral); i++) {
640 if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
641 gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
642 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
643 gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
645 for(i=0; (i<nrfp); i++)
647 push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line,wi);
651 void push_dihedraltype(directive d,t_params bt[],
652 t_bond_atomtype bat,char *line,
655 const char *formal[MAXATOMLIST+1] = {
664 const char *formnl[MAXATOMLIST+1] = {
670 "%*s%*s%*s%*s%*s%*s",
671 "%*s%*s%*s%*s%*s%*s%*s"
673 const char *formlf[MAXFORCEPARAM] = {
679 "%lf%lf%lf%lf%lf%lf",
680 "%lf%lf%lf%lf%lf%lf%lf",
681 "%lf%lf%lf%lf%lf%lf%lf%lf",
682 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
683 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
684 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
685 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
687 int i,ft,ftype,nn,nrfp,nrfpA,nrfpB,nral;
689 char alc[MAXATOMLIST+1][20];
690 double c[MAXFORCEPARAM];
692 gmx_bool bAllowRepeat;
695 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
697 * We first check for 2 atoms with the 3th column being an integer
698 * defining the type. If this isn't the case, we try it with 4 atoms
699 * and the 5th column defining the dihedral type.
701 nn=sscanf(line,formal[4],alc[0],alc[1],alc[2],alc[3],alc[4]);
702 if(nn>=3 && strlen(alc[2])==1 && isdigit(alc[2][0])) {
704 ft = strtol(alc[nral],NULL,10);
705 /* Move atom types around a bit and use 'X' for wildcard atoms
706 * to create a 4-atom dihedral definition with arbitrary atoms in
710 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
711 strcpy(alc[3],alc[1]);
714 /* alc[0] stays put */
716 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
718 strcpy(alc[2],alc[1]);
719 strcpy(alc[1],alc[0]);
722 } else if(nn==5 && strlen(alc[4])==1 && isdigit(alc[4][0])) {
724 ft = strtol(alc[nral],NULL,10);
726 sprintf(errbuf,"Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)",nn-1);
727 warning_error(wi,errbuf);
733 /* Previously, we have always overwritten parameters if e.g. a torsion
734 with the same atomtypes occurs on multiple lines. However, CHARMM and
735 some other force fields specify multiple dihedrals over some bonds,
736 including cosines with multiplicity 6 and somethimes even higher.
737 Thus, they cannot be represented with Ryckaert-Bellemans terms.
738 To add support for these force fields, Dihedral type 9 is identical to
739 normal proper dihedrals, but repeated entries are allowed.
746 bAllowRepeat = FALSE;
750 ftype = ifunc_index(d,ft);
752 nrfpA = interaction_function[ftype].nrfpA;
753 nrfpB = interaction_function[ftype].nrfpB;
755 strcpy(f1,formnl[nral]);
756 strcat(f1,formlf[nrfp-1]);
758 /* Check number of parameters given */
759 if ((nn=sscanf(line,f1,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5],&c[6],&c[7],&c[8],&c[9],&c[10],&c[11]))
762 /* Copy the B-state from the A-state */
763 copy_B_from_A(ftype,c);
766 warning_error(wi,"Not enough parameters");
767 } else if (nn > nrfpA && nn < nrfp) {
768 warning_error(wi,"Too many parameters or not enough parameters for topology B");
769 } else if (nn > nrfp) {
770 warning_error(wi,"Too many parameters");
772 for(i=nn; (i<nrfp); i++)
777 for(i=0; (i<4); i++) {
778 if(!strcmp(alc[i],"X"))
781 if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
782 gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
785 for(i=0; (i<nrfp); i++)
787 /* Always use 4 atoms here, since we created two wildcard atoms
788 * if there wasn't of them 4 already.
790 push_bondtype (&(bt[ftype]),&p,4,ftype,bAllowRepeat,line,wi);
794 void push_nbt(directive d,t_nbparam **nbt,gpp_atomtype_t atype,
795 char *pline,int nb_funct,
799 const char *form2="%*s%*s%*s%lf%lf";
800 const char *form3="%*s%*s%*s%lf%lf%lf";
801 const char *form4="%*s%*s%*s%lf%lf%lf%lf";
802 const char *form5="%*s%*s%*s%lf%lf%lf%lf%lf";
804 int i,f,n,ftype,atnr,nrfp;
812 if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
817 ftype=ifunc_index(d,f);
819 if (ftype != nb_funct) {
820 sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
821 interaction_function[ftype].longname,
822 interaction_function[nb_funct].longname);
823 warning_error(wi,errbuf);
827 /* Get the force parameters */
829 if (ftype == F_LJ14) {
830 n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
835 /* When the B topology parameters are not set,
836 * copy them from topology A
838 for(i=n; i<nrfp; i++)
841 else if (ftype == F_LJC14_Q) {
842 n = sscanf(pline,form5,&c[0],&c[1],&c[2],&c[3],&dum);
844 incorrect_n_param(wi);
848 else if (nrfp == 2) {
849 if (sscanf(pline,form3,&c[0],&c[1],&dum) != 2) {
850 incorrect_n_param(wi);
854 else if (nrfp == 3) {
855 if (sscanf(pline,form4,&c[0],&c[1],&c[2],&dum) != 3) {
856 incorrect_n_param(wi);
861 gmx_fatal(FARGS,"Number of force parameters for nonbonded interactions is %d"
862 " in file %s, line %d",nrfp,__FILE__,__LINE__);
864 for(i=0; (i<nrfp); i++)
867 /* Put the parameters in the matrix */
868 if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
869 gmx_fatal(FARGS,"Atomtype %s not found",a0);
870 if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
871 gmx_fatal(FARGS,"Atomtype %s not found",a1);
872 nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
876 for (i=0; i<nrfp; i++)
877 bId = bId && (nbp->c[i] == cr[i]);
879 sprintf(errbuf,"Overriding non-bonded parameters,");
881 fprintf(stderr," old:");
882 for (i=0; i<nrfp; i++)
883 fprintf(stderr," %g",nbp->c[i]);
884 fprintf(stderr," new\n%s\n",pline);
888 for (i=0; i<nrfp; i++)
893 push_gb_params (gpp_atomtype_t at, char *line,
898 double radius,vol,surftens,gb_radius,S_hct;
899 char atypename[STRLEN];
902 if ( (nfield = sscanf(line,"%s%lf%lf%lf%lf%lf",atypename,&radius,&vol,&surftens,&gb_radius,&S_hct)) != 6)
904 sprintf(errbuf,"Too few gb parameters for type %s\n",atypename);
908 /* Search for atomtype */
909 atype = get_atomtype_type(atypename,at);
913 printf("Couldn't find topology match for atomtype %s\n",atypename);
917 set_atomtype_gbparam(at,atype,radius,vol,surftens,gb_radius,S_hct);
921 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
922 t_bond_atomtype bat, char *line,
925 const char *formal="%s%s%s%s%s%s%s%s";
927 int i,j,ft,ftype,nn,nrfp,nrfpA,nrfpB;
929 int nxcmap,nycmap,ncmap,read_cmap,sl,nct;
930 char s[20],alc[MAXATOMLIST+2][20];
932 gmx_bool bAllowRepeat;
935 /* Keep the compiler happy */
939 if((nn=sscanf(line,formal,alc[0],alc[1],alc[2],alc[3],alc[4],alc[5],alc[6],alc[7])) != nral+3)
941 sprintf(errbuf,"Incorrect number of atomtypes for cmap (%d instead of 5)",nn-1);
942 warning_error(wi,errbuf);
946 /* Compute an offset for each line where the cmap parameters start
947 * ie. where the atom types and grid spacing information ends
951 start += (int)strlen(alc[i]);
954 /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
955 /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
956 start = start + nn -1;
958 ft = strtol(alc[nral],NULL,10);
959 nxcmap = strtol(alc[nral+1],NULL,10);
960 nycmap = strtol(alc[nral+2],NULL,10);
962 /* Check for equal grid spacing in x and y dims */
965 gmx_fatal(FARGS,"Not the same grid spacing in x and y for cmap grid: x=%d, y=%d",nxcmap,nycmap);
968 ncmap = nxcmap*nycmap;
969 ftype = ifunc_index(d,ft);
970 nrfpA = strtol(alc[6],NULL,10)*strtol(alc[6],NULL,10);
971 nrfpB = strtol(alc[7],NULL,10)*strtol(alc[7],NULL,10);
974 /* Allocate memory for the CMAP grid */
976 srenew(bt->cmap,bt->ncmap);
978 /* Read in CMAP parameters */
982 while(isspace(*(line+start+sl)))
986 nn=sscanf(line+start+sl," %s ",s);
988 bt->cmap[i+(bt->ncmap)-nrfp]=strtod(s,NULL);
996 gmx_fatal(FARGS,"Error in reading cmap parameter for angle %s %s %s %s %s",alc[0],alc[1],alc[2],alc[3],alc[4]);
1001 /* Check do that we got the number of parameters we expected */
1002 if(read_cmap==nrfpA)
1004 for(i=0;i<ncmap;i++)
1006 bt->cmap[i+ncmap]=bt->cmap[i];
1013 warning_error(wi,"Not enough cmap parameters");
1015 else if(read_cmap > nrfpA && read_cmap < nrfp)
1017 warning_error(wi,"Too many cmap parameters or not enough parameters for topology B");
1019 else if(read_cmap>nrfp)
1021 warning_error(wi,"Too many cmap parameters");
1026 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1027 * so we can safely assign them each time
1029 bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1030 bt->nc = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1031 nct = (nral+1) * bt->nc;
1033 /* Allocate memory for the cmap_types information */
1034 srenew(bt->cmap_types,nct);
1036 for(i=0; (i<nral); i++)
1038 if(at && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1040 gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
1042 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1044 gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
1047 /* Assign a grid number to each cmap_type */
1048 bt->cmap_types[bt->nct++]=get_bond_atomtype_type(alc[i],bat);
1051 /* Assign a type number to this cmap */
1052 bt->cmap_types[bt->nct++]=bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1054 /* Check for the correct number of atoms (again) */
1057 gmx_fatal(FARGS,"Incorrect number of atom types (%d) in cmap type %d\n",nct,bt->nc);
1060 /* Is this correct?? */
1061 for(i=0;i<MAXFORCEPARAM;i++)
1066 /* Push the bond to the bondlist */
1067 push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line,wi);
1071 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
1073 int type,char *ctype,int ptype,
1074 char *resnumberic,int cgnumber,
1075 char *resname,char *name,real m0,real q0,
1076 int typeB,char *ctypeB,real mB,real qB)
1078 int j,resind=0,resnr;
1082 if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1083 gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
1085 j = strlen(resnumberic) - 1;
1086 if (isdigit(resnumberic[j])) {
1089 ric = resnumberic[j];
1090 if (j == 0 || !isdigit(resnumberic[j-1])) {
1091 gmx_fatal(FARGS,"Invalid residue number '%s' for atom %d",
1092 resnumberic,atomnr);
1095 resnr = strtol(resnumberic,NULL,10);
1098 resind = at->atom[nr-1].resind;
1100 if (nr == 0 || strcmp(resname,*at->resinfo[resind].name) != 0 ||
1101 resnr != at->resinfo[resind].nr ||
1102 ric != at->resinfo[resind].ic) {
1108 at->nres = resind + 1;
1109 srenew(at->resinfo,at->nres);
1110 at->resinfo[resind].name = put_symtab(symtab,resname);
1111 at->resinfo[resind].nr = resnr;
1112 at->resinfo[resind].ic = ric;
1114 resind = at->atom[at->nr-1].resind;
1117 /* New atom instance
1118 * get new space for arrays
1120 srenew(at->atom,nr+1);
1121 srenew(at->atomname,nr+1);
1122 srenew(at->atomtype,nr+1);
1123 srenew(at->atomtypeB,nr+1);
1126 at->atom[nr].type = type;
1127 at->atom[nr].ptype = ptype;
1128 at->atom[nr].q = q0;
1129 at->atom[nr].m = m0;
1130 at->atom[nr].typeB = typeB;
1131 at->atom[nr].qB = qB;
1132 at->atom[nr].mB = mB;
1134 at->atom[nr].resind = resind;
1135 at->atom[nr].atomnumber = atomicnumber;
1136 at->atomname[nr] = put_symtab(symtab,name);
1137 at->atomtype[nr] = put_symtab(symtab,ctype);
1138 at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
1142 void push_cg(t_block *block, int *lastindex, int index, int a)
1145 fprintf (debug,"Index %d, Atom %d\n",index,a);
1147 if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
1148 /* add a new block */
1150 srenew(block->index,block->nr+1);
1152 block->index[block->nr] = a + 1;
1156 void push_atom(t_symtab *symtab,t_block *cgs,
1157 t_atoms *at,gpp_atomtype_t atype,char *line,int *lastcg,
1161 int cgnumber,atomnr,type,typeB,nscan;
1162 char id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
1163 resnumberic[STRLEN],resname[STRLEN],name[STRLEN],check[STRLEN];
1167 /* Make a shortcut for writing in this molecule */
1170 /* Fixed parameters */
1171 if (sscanf(line,"%s%s%s%s%s%d",
1172 id,ctype,resnumberic,resname,name,&cgnumber) != 6) {
1176 sscanf(id,"%d",&atomnr);
1177 if ((type = get_atomtype_type(ctype,atype)) == NOTSET)
1178 gmx_fatal(FARGS,"Atomtype %s not found",ctype);
1179 ptype = get_atomtype_ptype(type,atype);
1181 /* Set default from type */
1182 q0 = get_atomtype_qA(type,atype);
1183 m0 = get_atomtype_massA(type,atype);
1188 /* Optional parameters */
1189 nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1190 &q,&m,ctypeB,&qb,&mb,check);
1192 /* Nasty switch that falls thru all the way down! */
1198 if ((typeB = get_atomtype_type(ctypeB,atype)) == NOTSET) {
1199 gmx_fatal(FARGS,"Atomtype %s not found",ctypeB);
1201 qB = get_atomtype_qA(typeB,atype);
1202 mB = get_atomtype_massA(typeB,atype);
1208 warning_error(wi,"Too many parameters");
1215 fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
1217 push_cg(cgs,lastcg,cgnumber,nr);
1219 push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
1220 type,ctype,ptype,resnumberic,cgnumber,
1221 resname,name,m0,q0,typeB,
1222 typeB==type ? ctype : ctypeB,mB,qB);
1225 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line,
1232 if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
1233 warning_error(wi,"Expected a molecule type name and nrexcl");
1236 /* Test if this atomtype overwrites another */
1239 if (gmx_strcasecmp(*((*mol)[i].name),type) == 0)
1240 gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1246 newmol = &((*mol)[*nmol-1]);
1247 init_molinfo(newmol);
1249 /* Fill in the values */
1250 newmol->name = put_symtab(symtab,type);
1251 newmol->nrexcl = nrexcl;
1252 newmol->excl_set = FALSE;
1255 static gmx_bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1256 t_param *p,int c_start,gmx_bool bB,gmx_bool bGenPairs)
1258 int i,j,ti,tj,ntype;
1261 int nr = bt[ftype].nr;
1262 int nral = NRAL(ftype);
1263 int nrfp = interaction_function[ftype].nrfpA;
1264 int nrfpB= interaction_function[ftype].nrfpB;
1266 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1271 /* First test the generated-pair position to save
1272 * time when we have 1000*1000 entries for e.g. OPLS...
1276 ti=at->atom[p->a[0]].typeB;
1277 tj=at->atom[p->a[1]].typeB;
1279 ti=at->atom[p->a[0]].type;
1280 tj=at->atom[p->a[1]].type;
1282 pi=&(bt[ftype].param[ntype*ti+tj]);
1283 bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1286 /* Search explicitly if we didnt find it */
1288 for (i=0; ((i < nr) && !bFound); i++) {
1289 pi=&(bt[ftype].param[i]);
1291 for (j=0; ((j < nral) &&
1292 (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1294 for (j=0; ((j < nral) &&
1295 (at->atom[p->a[j]].type == pi->a[j])); j++);
1302 if (nrfp+nrfpB > MAXFORCEPARAM) {
1303 gmx_incons("Too many force parameters");
1305 for (j=c_start; (j < nrfpB); j++)
1306 p->c[nrfp+j] = pi->c[j];
1309 for (j=c_start; (j < nrfp); j++)
1313 for (j=c_start; (j < nrfp); j++)
1319 static gmx_bool default_cmap_params(int ftype, t_params bondtype[],
1320 t_atoms *at, gpp_atomtype_t atype,
1321 t_param *p, gmx_bool bB,
1322 int *cmap_type, int *nparam_def)
1324 int i,j,nparam_found;
1326 gmx_bool bFound=FALSE;
1331 /* Match the current cmap angle against the list of cmap_types */
1332 for(i=0;i<bondtype->nct && !bFound;i+=6)
1341 (get_atomtype_batype(at->atom[p->a[0]].type,atype)==bondtype->cmap_types[i]) &&
1342 (get_atomtype_batype(at->atom[p->a[1]].type,atype)==bondtype->cmap_types[i+1]) &&
1343 (get_atomtype_batype(at->atom[p->a[2]].type,atype)==bondtype->cmap_types[i+2]) &&
1344 (get_atomtype_batype(at->atom[p->a[3]].type,atype)==bondtype->cmap_types[i+3]) &&
1345 (get_atomtype_batype(at->atom[p->a[4]].type,atype)==bondtype->cmap_types[i+4]))
1347 /* Found cmap torsion */
1349 ct = bondtype->cmap_types[i+5];
1355 /* If we did not find a matching type for this cmap torsion */
1358 gmx_fatal(FARGS,"Unknown cmap torsion between atoms %d %d %d %d %d\n",
1359 p->a[0]+1,p->a[1]+1,p->a[2]+1,p->a[3]+1,p->a[4]+1);
1362 *nparam_def = nparam_found;
1368 static gmx_bool default_params(int ftype,t_params bt[],
1369 t_atoms *at,gpp_atomtype_t atype,
1370 t_param *p,gmx_bool bB,
1371 t_param **param_def,
1374 int i,j,nparam_found;
1375 gmx_bool bFound,bSame;
1378 int nr = bt[ftype].nr;
1379 int nral = NRAL(ftype);
1380 int nrfpA= interaction_function[ftype].nrfpA;
1381 int nrfpB= interaction_function[ftype].nrfpB;
1383 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1387 /* We allow wildcards now. The first type (with or without wildcards) that
1388 * fits is used, so you should probably put the wildcarded bondtypes
1389 * at the end of each section.
1393 /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1394 * special case for this. Check for B state outside loop to speed it up.
1396 if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1400 for (i=0; ((i < nr) && !bFound); i++)
1402 pi=&(bt[ftype].param[i]);
1405 ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1406 ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1407 ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1408 ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1415 for (i=0; ((i < nr) && !bFound); i++)
1417 pi=&(bt[ftype].param[i]);
1420 ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&
1421 ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1422 ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1423 ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1427 /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1428 * The rules in that case is that additional matches HAVE to be on adjacent lines!
1434 /* Continue from current i value */
1435 for(j=i+1 ; j<nr && bSame ; j+=2)
1437 pj=&(bt[ftype].param[j]);
1438 bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1441 /* nparam_found will be increased as long as the numbers match */
1444 } else { /* Not a dihedral */
1445 for (i=0; ((i < nr) && !bFound); i++) {
1446 pi=&(bt[ftype].param[i]);
1448 for (j=0; ((j < nral) &&
1449 (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1452 for (j=0; ((j < nral) &&
1453 (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1462 *nparam_def = nparam_found;
1469 void push_bond(directive d,t_params bondtype[],t_params bond[],
1470 t_atoms *at,gpp_atomtype_t atype,char *line,
1471 gmx_bool bBonded,gmx_bool bGenPairs,real fudgeQQ,
1472 gmx_bool bZero,gmx_bool *bWarn_copy_A_B,
1475 const char *aaformat[MAXATOMLIST]= {
1483 const char *asformat[MAXATOMLIST]= {
1488 "%*s%*s%*s%*s%*s%*s",
1489 "%*s%*s%*s%*s%*s%*s%*s"
1491 const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1492 int nr,i,j,nral,nral_fmt,nread,ftype;
1493 char format[STRLEN];
1494 /* One force parameter more, so we can check if we read too many */
1495 double cc[MAXFORCEPARAM+1];
1496 int aa[MAXATOMLIST+1];
1497 t_param param,paramB,*param_defA,*param_defB;
1498 gmx_bool bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1499 int nparam_defA,nparam_defB;
1502 nparam_defA=nparam_defB=0;
1504 ftype = ifunc_index(d,1);
1506 for(j=0; j<MAXATOMLIST; j++)
1510 bDef = (NRFP(ftype) > 0);
1512 if (ftype == F_SETTLE)
1514 /* SETTLE acts on 3 atoms, but the topology format only specifies
1515 * the first atom (for historical reasons).
1524 nread = sscanf(line,aaformat[nral_fmt-1],
1525 &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1527 if (ftype == F_SETTLE)
1534 if (nread < nral_fmt)
1539 else if (nread > nral_fmt)
1541 /* this is a hack to allow for virtual sites with swapped parity */
1542 bSwapParity = (aa[nral]<0);
1545 aa[nral] = -aa[nral];
1547 ftype = ifunc_index(d,aa[nral]);
1556 gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1557 interaction_function[F_VSITE3FAD].longname,
1558 interaction_function[F_VSITE3OUT].longname);
1564 /* Check for double atoms and atoms out of bounds */
1565 for(i=0; (i<nral); i++) {
1566 if ( aa[i] < 1 || aa[i] > at->nr )
1567 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1568 "Atom index (%d) in %s out of bounds (1-%d).\n"
1569 "This probably means that you have inserted topology section \"%s\"\n"
1570 "in a part belonging to a different molecule than you intended to.\n"
1571 "In that case move the \"%s\" section to the right molecule.",
1572 get_warning_file(wi),get_warning_line(wi),
1573 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1574 for(j=i+1; (j<nral); j++)
1575 if (aa[i] == aa[j]) {
1576 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1581 /* default force parameters */
1582 for(j=0; (j<MAXATOMLIST); j++)
1583 param.a[j] = aa[j]-1;
1584 for(j=0; (j<MAXFORCEPARAM); j++)
1587 /* Get force params for normal and free energy perturbation
1588 * studies, as determined by types!
1592 bFoundA = default_params(ftype,bondtype,at,atype,¶m,FALSE,¶m_defA,&nparam_defA);
1594 /* Copy the A-state and B-state default parameters */
1595 for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1596 param.c[j] = param_defA->c[j];
1598 bFoundB = default_params(ftype,bondtype,at,atype,¶m,TRUE,¶m_defB,&nparam_defB);
1600 /* Copy only the B-state default parameters */
1601 for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1602 param.c[j] = param_defB->c[j];
1604 } else if (ftype == F_LJ14) {
1605 bFoundA = default_nb_params(ftype, bondtype,at,¶m,0,FALSE,bGenPairs);
1606 bFoundB = default_nb_params(ftype, bondtype,at,¶m,0,TRUE, bGenPairs);
1607 } else if (ftype == F_LJC14_Q) {
1608 param.c[0] = fudgeQQ;
1609 /* Fill in the A-state charges as default parameters */
1610 param.c[1] = at->atom[param.a[0]].q;
1611 param.c[2] = at->atom[param.a[1]].q;
1612 /* The default LJ parameters are the standard 1-4 parameters */
1613 bFoundA = default_nb_params(F_LJ14,bondtype,at,¶m,3,FALSE,bGenPairs);
1615 } else if (ftype == F_LJC_PAIRS_NB) {
1616 /* Defaults are not supported here */
1620 gmx_incons("Unknown function type in push_bond");
1623 if (nread > nral_fmt) {
1624 /* Manually specified parameters - in this case we discard multiple torsion info! */
1626 strcpy(format,asformat[nral_fmt-1]);
1627 strcat(format,ccformat);
1629 nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1630 &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1632 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1634 /* We only have to issue a warning if these atoms are perturbed! */
1636 for(j=0; (j<nral); j++)
1637 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1639 if (bPert && *bWarn_copy_A_B)
1642 "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1644 *bWarn_copy_A_B = FALSE;
1647 /* If only the A parameters were specified, copy them to the B state */
1648 /* The B-state parameters correspond to the first nrfpB
1649 * A-state parameters.
1651 for(j=0; (j<NRFPB(ftype)); j++)
1652 cc[nread++] = cc[j];
1655 /* If nread was 0 or EOF, no parameters were read => use defaults.
1656 * If nread was nrfpA we copied above so nread=nrfp.
1657 * If nread was nrfp we are cool.
1658 * For F_LJC14_Q we allow supplying fudgeQQ only.
1659 * Anything else is an error!
1661 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1662 !(ftype == F_LJC14_Q && nread == 1))
1664 gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1665 nread,NRFPA(ftype),NRFP(ftype),
1666 interaction_function[ftype].longname);
1669 for(j=0; (j<nread); j++)
1672 /* Check whether we have to use the defaults */
1673 if (nread == NRFP(ftype))
1680 /* nread now holds the number of force parameters read! */
1685 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1688 if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1691 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1692 "Please specify perturbed parameters manually for this torsion in your topology!");
1693 warning_error(wi,errbuf);
1697 if (nread > 0 && nread < NRFPA(ftype))
1699 /* Issue an error, do not use defaults */
1700 sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1701 warning_error(wi,errbuf);
1704 if (nread == 0 || nread == EOF)
1708 if (interaction_function[ftype].flags & IF_VSITE)
1710 /* set them to NOTSET, will be calculated later */
1711 for(j=0; (j<MAXFORCEPARAM); j++)
1712 param.c[j] = NOTSET;
1715 param.C1 = -1; /* flag to swap parity of vsite construction */
1721 fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1722 interaction_function[ftype].longname);
1726 sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1727 warning_error(wi,errbuf);
1738 param.C0 = 360 - param.C0;
1741 param.C2 = -param.C2;
1748 /* We only have to issue a warning if these atoms are perturbed! */
1750 for(j=0; (j<nral); j++)
1751 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1755 sprintf(errbuf,"No default %s types for perturbed atoms, "
1756 "using normal values",interaction_function[ftype].longname);
1763 if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1764 && param.c[5]!=param.c[2])
1765 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1766 " %s multiplicity can not be perturbed %f!=%f",
1767 get_warning_file(wi),get_warning_line(wi),
1768 interaction_function[ftype].longname,
1769 param.c[2],param.c[5]);
1771 if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1772 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1773 " %s table number can not be perturbed %d!=%d",
1774 get_warning_file(wi),get_warning_line(wi),
1775 interaction_function[ftype].longname,
1776 (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1778 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1779 if (ftype==F_RBDIHS) {
1781 for(i=0;i<NRFP(ftype);i++) {
1789 /* Put the values in the appropriate arrays */
1790 add_param_to_list (&bond[ftype],¶m);
1792 /* Push additional torsions from FF for ftype==9 if we have them.
1793 * We have already checked that the A/B states do not differ in this case,
1794 * so we do not have to double-check that again, or the vsite stuff.
1795 * In addition, those torsions cannot be automatically perturbed.
1797 if(bDef && ftype==F_PDIHS)
1799 for(i=1;i<nparam_defA;i++)
1801 /* Advance pointer! */
1803 for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1804 param.c[j] = param_defA->c[j];
1805 /* And push the next term for this torsion */
1806 add_param_to_list (&bond[ftype],¶m);
1811 void push_cmap(directive d, t_params bondtype[], t_params bond[],
1812 t_atoms *at, gpp_atomtype_t atype, char *line,
1813 gmx_bool *bWarn_copy_A_B,
1816 const char *aaformat[MAXATOMLIST+1]=
1827 int i,j,nr,ftype,nral,nread,ncmap_params;
1829 int aa[MAXATOMLIST+1];
1832 t_param param,paramB,*param_defA,*param_defB;
1834 ftype = ifunc_index(d,1);
1836 nr = bondtype[ftype].nr;
1839 nread = sscanf(line,aaformat[nral-1],
1840 &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1847 else if (nread == nral)
1849 ftype = ifunc_index(d,1);
1852 /* Check for double atoms and atoms out of bounds */
1855 if(aa[i] < 1 || aa[i] > at->nr)
1857 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1858 "Atom index (%d) in %s out of bounds (1-%d).\n"
1859 "This probably means that you have inserted topology section \"%s\"\n"
1860 "in a part belonging to a different molecule than you intended to.\n"
1861 "In that case move the \"%s\" section to the right molecule.",
1862 get_warning_file(wi),get_warning_line(wi),
1863 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1866 for(j=i+1; (j<nral); j++)
1870 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1876 /* default force parameters */
1877 for(j=0; (j<MAXATOMLIST); j++)
1879 param.a[j] = aa[j]-1;
1881 for(j=0; (j<MAXFORCEPARAM); j++)
1886 /* Get the cmap type for this cmap angle */
1887 bFound = default_cmap_params(ftype,bondtype,at,atype,¶m,FALSE,&cmap_type,&ncmap_params);
1889 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
1890 if(bFound && ncmap_params==1)
1892 /* Put the values in the appropriate arrays */
1893 param.c[0]=cmap_type;
1894 add_param_to_list(&bond[ftype],¶m);
1898 /* This is essentially the same check as in default_cmap_params() done one more time */
1899 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
1900 param.a[0]+1,param.a[1]+1,param.a[2]+1,param.a[3]+1,param.a[4]+1);
1906 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1907 t_atoms *at,gpp_atomtype_t atype,char *line,
1911 int type,ftype,j,n,ret,nj,a;
1913 double *weight=NULL,weight_tot;
1916 /* default force parameters */
1917 for(j=0; (j<MAXATOMLIST); j++)
1918 param.a[j] = NOTSET;
1919 for(j=0; (j<MAXFORCEPARAM); j++)
1923 ret = sscanf(ptr,"%d%n",&a,&n);
1926 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1927 " Expected an atom index in section \"%s\"",
1928 get_warning_file(wi),get_warning_line(wi),
1933 ret = sscanf(ptr,"%d%n",&type,&n);
1935 ftype = ifunc_index(d,type);
1940 ret = sscanf(ptr,"%d%n",&a,&n);
1945 srenew(weight,nj+20);
1953 /* Here we use the A-state mass as a parameter.
1954 * Note that the B-state mass has no influence.
1956 weight[nj] = at->atom[atc[nj]].m;
1960 ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1963 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1964 " No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1965 get_warning_file(wi),get_warning_line(wi),
1969 gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1971 weight_tot += weight[nj];
1977 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1978 " Expected more than one atom index in section \"%s\"",
1979 get_warning_file(wi),get_warning_line(wi),
1982 if (weight_tot == 0)
1983 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1984 " The total mass of the construting atoms is zero",
1985 get_warning_file(wi),get_warning_line(wi));
1987 for(j=0; j<nj; j++) {
1988 param.a[1] = atc[j];
1990 param.c[1] = weight[j]/weight_tot;
1991 /* Put the values in the appropriate arrays */
1992 add_param_to_list (&bond[ftype],¶m);
1999 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
2007 if (sscanf(pline,"%s%d",type,&copies) != 2) {
2012 /* search moleculename */
2013 for (i=0; ((i<nrmols) && gmx_strcasecmp(type,*(mols[i].name))); i++)
2020 gmx_fatal(FARGS,"No such moleculetype %s",type);
2023 void init_block2(t_block2 *b2, int natom)
2028 snew(b2->nra,b2->nr);
2030 for(i=0; (i<b2->nr); i++)
2034 void done_block2(t_block2 *b2)
2039 for(i=0; (i<b2->nr); i++)
2047 void push_excl(char *line, t_block2 *b2)
2051 char base[STRLEN],format[STRLEN];
2053 if (sscanf(line,"%d",&i)==0)
2056 if ((1 <= i) && (i <= b2->nr))
2060 fprintf(debug,"Unbound atom %d\n",i-1);
2065 strcpy(format,base);
2066 strcat(format,"%d");
2067 n=sscanf(line,format,&j);
2069 if ((1 <= j) && (j <= b2->nr)) {
2071 srenew(b2->a[i],++(b2->nra[i]));
2072 b2->a[i][b2->nra[i]-1]=j;
2073 /* also add the reverse exclusion! */
2074 srenew(b2->a[j],++(b2->nra[j]));
2075 b2->a[j][b2->nra[j]-1]=i;
2079 gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
2084 void b_to_b2(t_blocka *b, t_block2 *b2)
2089 for(i=0; (i<b->nr); i++)
2090 for(j=b->index[i]; (j<b->index[i+1]); j++) {
2092 srenew(b2->a[i],++b2->nra[i]);
2093 b2->a[i][b2->nra[i]-1]=a;
2097 void b2_to_b(t_block2 *b2, t_blocka *b)
2103 for(i=0; (i<b2->nr); i++) {
2105 for(j=0; (j<b2->nra[i]); j++)
2106 b->a[nra+j]=b2->a[i][j];
2109 /* terminate list */
2113 static int icomp(const void *v1, const void *v2)
2115 return (*((atom_id *) v1))-(*((atom_id *) v2));
2118 void merge_excl(t_blocka *excl, t_block2 *b2)
2126 else if (b2->nr != excl->nr) {
2127 gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2131 fprintf(debug,"Entering merge_excl\n");
2133 /* First copy all entries from excl to b2 */
2136 /* Count and sort the exclusions */
2138 for(i=0; (i<b2->nr); i++) {
2139 if (b2->nra[i] > 0) {
2140 /* remove double entries */
2141 qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
2143 for(j=1; (j<b2->nra[i]); j++)
2144 if (b2->a[i][j]!=b2->a[i][k-1]) {
2145 b2->a[i][k]=b2->a[i][j];
2153 srenew(excl->a,excl->nra);
2158 int add_atomtype_decoupled(t_symtab *symtab,gpp_atomtype_t at,
2159 t_nbparam ***nbparam,t_nbparam ***pair)
2165 /* Define an atom type with all parameters set to zero (no interactions) */
2168 /* Type for decoupled atoms could be anything,
2169 * this should be changed automatically later when required.
2171 atom.ptype = eptAtom;
2172 for (i=0; (i<MAXFORCEPARAM); i++)
2175 nr = add_atomtype(at,symtab,&atom,"decoupled",¶m,-1,0.0,0.0,0.0,0,0,0);
2177 /* Add space in the non-bonded parameters matrix */
2178 realloc_nb_params(at,nbparam,pair);
2183 static void convert_pairs_to_pairsQ(t_params *plist,
2184 real fudgeQQ,t_atoms *atoms)
2190 /* Copy the pair list to the pairQ list */
2191 plist[F_LJC14_Q] = plist[F_LJ14];
2192 /* Empty the pair list */
2193 plist[F_LJ14].nr = 0;
2194 plist[F_LJ14].param = NULL;
2195 param = plist[F_LJC14_Q].param;
2196 for(i=0; i<plist[F_LJC14_Q].nr; i++) {
2199 param[i].c[0] = fudgeQQ;
2200 param[i].c[1] = atoms->atom[param[i].a[0]].q;
2201 param[i].c[2] = atoms->atom[param[i].a[1]].q;
2207 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
2216 atom = mol->atoms.atom;
2218 ntype = sqrt(nbp->nr);
2220 for (i=0; i<MAXATOMLIST; i++)
2221 param.a[i] = NOTSET;
2222 for (i=0; i<MAXFORCEPARAM; i++)
2223 param.c[i] = NOTSET;
2225 /* Add a pair interaction for all non-excluded atom pairs */
2227 for(i=0; i<n; i++) {
2228 for(j=i+1; j<n; j++) {
2230 for(k=excl->index[i]; k<excl->index[i+1]; k++) {
2231 if (excl->a[k] == j)
2235 if (nb_funct != F_LJ)
2236 gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2239 param.c[0] = atom[i].q;
2240 param.c[1] = atom[j].q;
2241 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2242 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2243 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],¶m);
2249 static void set_excl_all(t_blocka *excl)
2253 /* Get rid of the current exclusions and exclude all atom pairs */
2255 excl->nra = nat*nat;
2256 srenew(excl->a,excl->nra);
2258 for(i=0; i<nat; i++) {
2260 for(j=0; j<nat; j++)
2263 excl->index[nat] = k;
2266 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
2267 int couple_lam0,int couple_lam1)
2271 for(i=0; i<atoms->nr; i++) {
2272 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) {
2273 atoms->atom[i].q = 0.0;
2275 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) {
2276 atoms->atom[i].type = atomtype_decouple;
2278 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) {
2279 atoms->atom[i].qB = 0.0;
2281 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) {
2282 atoms->atom[i].typeB = atomtype_decouple;
2287 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
2288 int couple_lam0,int couple_lam1,
2289 gmx_bool bCoupleIntra,int nb_funct,t_params *nbp)
2291 convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
2292 if (!bCoupleIntra) {
2293 generate_LJCpairsNB(mol,nb_funct,nbp);
2294 set_excl_all(&mol->excls);
2296 decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);