1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
53 #include "gmx_fatal.h"
55 #include "gpp_atomtype.h"
56 #include "gpp_bond_atomtype.h"
58 void generate_nbparams(int comb,int ftype,t_params *plist,gpp_atomtype_t atype,
63 real c,bi,bj,ci,cj,ci0,ci1,ci2,cj0,cj1,cj2;
66 /* Lean mean shortcuts */
67 nr = get_atomtype_ntypes(atype);
69 snew(plist->param,nr*nr);
72 /* Fill the matrix with force parameters */
78 for(i=k=0; (i<nr); i++) {
79 for(j=0; (j<nr); j++,k++) {
80 for(nf=0; (nf<nrfp); nf++) {
81 ci = get_atomtype_nbparam(i,nf,atype);
82 cj = get_atomtype_nbparam(j,nf,atype);
84 plist->param[k].c[nf] = c;
90 case eCOMB_ARITHMETIC:
91 /* c0 and c1 are epsilon and sigma */
92 for (i=k=0; (i < nr); i++)
93 for (j=0; (j < nr); j++,k++) {
94 ci0 = get_atomtype_nbparam(i,0,atype);
95 cj0 = get_atomtype_nbparam(j,0,atype);
96 ci1 = get_atomtype_nbparam(i,1,atype);
97 cj1 = get_atomtype_nbparam(j,1,atype);
98 plist->param[k].c[0] = (ci0+cj0)*0.5;
99 plist->param[k].c[1] = sqrt(ci1*cj1);
103 case eCOMB_GEOM_SIG_EPS:
104 /* c0 and c1 are epsilon and sigma */
105 for (i=k=0; (i < nr); i++)
106 for (j=0; (j < nr); j++,k++) {
107 ci0 = get_atomtype_nbparam(i,0,atype);
108 cj0 = get_atomtype_nbparam(j,0,atype);
109 ci1 = get_atomtype_nbparam(i,1,atype);
110 cj1 = get_atomtype_nbparam(j,1,atype);
111 plist->param[k].c[0] = sqrt(ci0*cj0);
112 plist->param[k].c[1] = sqrt(ci1*cj1);
117 gmx_fatal(FARGS,"No such combination rule %d",comb);
120 gmx_incons("Topology processing, generate nb parameters");
124 /* Buckingham rules */
125 for(i=k=0; (i<nr); i++)
126 for(j=0; (j<nr); j++,k++) {
127 ci0 = get_atomtype_nbparam(i,0,atype);
128 cj0 = get_atomtype_nbparam(j,0,atype);
129 ci2 = get_atomtype_nbparam(i,2,atype);
130 cj2 = get_atomtype_nbparam(j,2,atype);
131 bi = get_atomtype_nbparam(i,1,atype);
132 bj = get_atomtype_nbparam(j,1,atype);
133 plist->param[k].c[0] = sqrt(ci0 * cj0);
134 if ((bi == 0) || (bj == 0))
135 plist->param[k].c[1] = 0;
137 plist->param[k].c[1] = 2.0/(1/bi+1/bj);
138 plist->param[k].c[2] = sqrt(ci2 * cj2);
143 sprintf(errbuf,"Invalid nonbonded type %s",
144 interaction_function[ftype].longname);
145 warning_error(wi,errbuf);
149 static void realloc_nb_params(gpp_atomtype_t at,
150 t_nbparam ***nbparam,t_nbparam ***pair)
152 /* Add space in the non-bonded parameters matrix */
153 int atnr = get_atomtype_ntypes(at);
154 srenew(*nbparam,atnr);
155 snew((*nbparam)[atnr-1],atnr);
158 snew((*pair)[atnr-1],atnr);
162 static void copy_B_from_A(int ftype,double *c)
166 nrfpA = NRFPA(ftype);
167 nrfpB = NRFPB(ftype);
169 /* Copy the B parameters from the first nrfpB A parameters */
170 for(i=0; (i<nrfpB); i++) {
175 void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
176 char *line,int nb_funct,
177 t_nbparam ***nbparam,t_nbparam ***pair,
184 t_xlate xl[eptNR] = {
192 int nr,i,nfields,j,pt,nfp0=-1;
194 char type[STRLEN],btype[STRLEN],ptype[STRLEN];
196 double c[MAXFORCEPARAM];
197 double radius,vol,surftens,gb_radius,S_hct;
198 char tmpfield[12][100]; /* Max 12 fields of width 100 */
203 gmx_bool have_atomic_number;
204 gmx_bool have_bonded_type;
209 /* First assign input line to temporary array */
210 nfields=sscanf(line,"%s%s%s%s%s%s%s%s%s%s%s%s",
211 tmpfield[0],tmpfield[1],tmpfield[2],tmpfield[3],tmpfield[4],tmpfield[5],
212 tmpfield[6],tmpfield[7],tmpfield[8],tmpfield[9],tmpfield[10],tmpfield[11]);
214 /* Comments on optional fields in the atomtypes section:
216 * The force field format is getting a bit old. For OPLS-AA we needed
217 * to add a special bonded atomtype, and for Gerrit Groenhofs QM/MM stuff
218 * we also needed the atomic numbers.
219 * To avoid making all old or user-generated force fields unusable we
220 * have introduced both these quantities as optional columns, and do some
221 * acrobatics to check whether they are present or not.
222 * This will all look much nicer when we switch to XML... sigh.
224 * Field 0 (mandatory) is the nonbonded type name. (string)
225 * Field 1 (optional) is the bonded type (string)
226 * Field 2 (optional) is the atomic number (int)
227 * Field 3 (mandatory) is the mass (numerical)
228 * Field 4 (mandatory) is the charge (numerical)
229 * Field 5 (mandatory) is the particle type (single character)
230 * This is followed by a number of nonbonded parameters.
232 * The safest way to identify the format is the particle type field.
234 * So, here is what we do:
236 * A. Read in the first six fields as strings
237 * B. If field 3 (starting from 0) is a single char, we have neither
238 * bonded_type or atomic numbers.
239 * C. If field 5 is a single char we have both.
240 * D. If field 4 is a single char we check field 1. If this begins with
241 * an alphabetical character we have bonded types, otherwise atomic numbers.
250 if( (strlen(tmpfield[5])==1) && isalpha(tmpfield[5][0]) )
252 have_bonded_type = TRUE;
253 have_atomic_number = TRUE;
255 else if( (strlen(tmpfield[3])==1) && isalpha(tmpfield[3][0]) )
257 have_bonded_type = FALSE;
258 have_atomic_number = FALSE;
262 have_bonded_type = ( isalpha(tmpfield[1][0]) != 0 );
263 have_atomic_number = !have_bonded_type;
266 /* optional fields */
279 if( have_atomic_number )
281 if ( have_bonded_type )
283 nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
284 type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],
285 &radius,&vol,&surftens,&gb_radius);
294 /* have_atomic_number && !have_bonded_type */
295 nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf",
296 type,&atomnr,&m,&q,ptype,&c[0],&c[1],
297 &radius,&vol,&surftens,&gb_radius);
307 if ( have_bonded_type )
309 /* !have_atomic_number && have_bonded_type */
310 nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
311 type,btype,&m,&q,ptype,&c[0],&c[1],
312 &radius,&vol,&surftens,&gb_radius);
321 /* !have_atomic_number && !have_bonded_type */
322 nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf",
323 type,&m,&q,ptype,&c[0],&c[1],
324 &radius,&vol,&surftens,&gb_radius);
333 if( !have_bonded_type )
338 if( !have_atomic_number )
348 if( have_atomic_number )
350 if ( have_bonded_type )
352 nread = sscanf(line,"%s%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
353 type,btype,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
354 &radius,&vol,&surftens,&gb_radius);
363 /* have_atomic_number && !have_bonded_type */
364 nread = sscanf(line,"%s%d%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
365 type,&atomnr,&m,&q,ptype,&c[0],&c[1],&c[2],
366 &radius,&vol,&surftens,&gb_radius);
376 if ( have_bonded_type )
378 /* !have_atomic_number && have_bonded_type */
379 nread = sscanf(line,"%s%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
380 type,btype,&m,&q,ptype,&c[0],&c[1],&c[2],
381 &radius,&vol,&surftens,&gb_radius);
390 /* !have_atomic_number && !have_bonded_type */
391 nread = sscanf(line,"%s%lf%lf%s%lf%lf%lf%lf%lf%lf%lf",
392 type,&m,&q,ptype,&c[0],&c[1],&c[2],
393 &radius,&vol,&surftens,&gb_radius);
402 if( !have_bonded_type )
407 if( !have_atomic_number )
415 gmx_fatal(FARGS,"Invalid function type %d in push_at %s %d",nb_funct,
418 for(j=nfp0; (j<MAXFORCEPARAM); j++)
421 if(strlen(type)==1 && isdigit(type[0]))
422 gmx_fatal(FARGS,"Atom type names can't be single digits.");
424 if(strlen(btype)==1 && isdigit(btype[0]))
425 gmx_fatal(FARGS,"Bond atom type names can't be single digits.");
427 /* Hack to read old topologies */
428 if (gmx_strcasecmp(ptype,"D") == 0)
430 for(j=0; (j<eptNR); j++)
431 if (gmx_strcasecmp(ptype,xl[j].entry) == 0)
434 gmx_fatal(FARGS,"Invalid particle type %s on line %s",
438 fprintf(debug,"ptype: %s\n",ptype_str[pt]);
443 for (i=0; (i<MAXFORCEPARAM); i++)
446 if ((batype_nr = get_bond_atomtype_type(btype,bat)) == NOTSET)
447 add_bond_atomtype(bat,symtab,btype);
448 batype_nr = get_bond_atomtype_type(btype,bat);
450 if ((nr = get_atomtype_type(type,at)) != NOTSET) {
451 sprintf(errbuf,"Overriding atomtype %s",type);
453 if ((nr = set_atomtype(nr,at,symtab,atom,type,param,batype_nr,
454 radius,vol,surftens,atomnr,gb_radius,S_hct)) == NOTSET)
455 gmx_fatal(FARGS,"Replacing atomtype %s failed",type);
457 else if ((nr = add_atomtype(at,symtab,atom,type,param,
458 batype_nr,radius,vol,
459 surftens,atomnr,gb_radius,S_hct)) == NOTSET)
460 gmx_fatal(FARGS,"Adding atomtype %s failed",type);
462 /* Add space in the non-bonded parameters matrix */
463 realloc_nb_params(at,nbparam,pair);
469 static void push_bondtype(t_params * bt,
473 gmx_bool bAllowRepeat,
478 gmx_bool bTest,bFound,bCont,bId;
480 int nrfp = NRFP(ftype);
483 /* If bAllowRepeat is TRUE, we allow multiple entries as long as they
484 are on directly _adjacent_ lines.
487 /* First check if our atomtypes are _identical_ (not reversed) to the previous
488 entry. If they are not identical we search for earlier duplicates. If they are
489 we can skip it, since we already searched for the first line
496 if(bAllowRepeat && nr>1)
498 for (j=0,bCont=TRUE; (j < nral); j++)
500 bCont=bCont && (b->a[j] == bt->param[nr-2].a[j]);
504 /* Search for earlier duplicates if this entry was not a continuation
505 from the previous line.
510 for (i=0; (i < nr); i++) {
512 for (j=0; (j < nral); j++)
513 bTest=(bTest && (b->a[j] == bt->param[i].a[j]));
516 for(j=0; (j<nral); j++)
517 bTest=(bTest && (b->a[nral-1-j] == bt->param[i].a[j]));
522 for (j=0; (j < nrfp); j++)
523 bId = bId && (bt->param[i].c[j] == b->c[j]);
525 sprintf(errbuf,"Overriding %s parameters.%s",
526 interaction_function[ftype].longname,
527 (ftype==F_PDIHS) ? "\nUse dihedraltype 4 to allow several multiplicity terms." : "");
529 fprintf(stderr," old:");
530 for (j=0; (j < nrfp); j++)
531 fprintf(stderr," %g",bt->param[i].c[j]);
532 fprintf(stderr," \n new: %s\n\n",line);
536 for (j=0; (j < nrfp); j++)
537 bt->param[i].c[j] = b->c[j];
546 /* fill the arrays up and down */
547 memcpy(bt->param[bt->nr].c, b->c,sizeof(b->c));
548 memcpy(bt->param[bt->nr].a, b->a,sizeof(b->a));
549 memcpy(bt->param[bt->nr+1].c,b->c,sizeof(b->c));
551 /* The definitions of linear angles depend on the order of atoms,
552 * that means that for atoms i-j-k, with certain parameter a, the
553 * corresponding k-j-i angle will have parameter 1-a.
555 if (ftype == F_LINEAR_ANGLES)
557 bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0];
558 bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2];
561 for (j=0; (j < nral); j++)
563 bt->param[bt->nr+1].a[j] = b->a[nral-1-j];
570 void push_bt(directive d,t_params bt[],int nral,
572 t_bond_atomtype bat,char *line,
575 const char *formal[MAXATOMLIST+1] = {
584 const char *formnl[MAXATOMLIST+1] = {
590 "%*s%*s%*s%*s%*s%*s",
591 "%*s%*s%*s%*s%*s%*s%*s"
593 const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
594 int i,ft,ftype,nn,nrfp,nrfpA,nrfpB;
596 char alc[MAXATOMLIST+1][20];
597 /* One force parameter more, so we can check if we read too many */
598 double c[MAXFORCEPARAM+1];
602 if ((bat && at) || (!bat && !at))
603 gmx_incons("You should pass either bat or at to push_bt");
605 /* Make format string (nral ints+functype) */
606 if ((nn=sscanf(line,formal[nral],
607 alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
608 sprintf(errbuf,"Not enough atomtypes (%d instead of %d)",nn-1,nral);
609 warning_error(wi,errbuf);
613 ft = strtol(alc[nral],NULL,10);
614 ftype = ifunc_index(d,ft);
616 nrfpA = interaction_function[ftype].nrfpA;
617 nrfpB = interaction_function[ftype].nrfpB;
618 strcpy(f1,formnl[nral]);
620 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]))
623 /* Copy the B-state from the A-state */
624 copy_B_from_A(ftype,c);
627 warning_error(wi,"Not enough parameters");
628 } else if (nn > nrfpA && nn < nrfp) {
629 warning_error(wi,"Too many parameters or not enough parameters for topology B");
630 } else if (nn > nrfp) {
631 warning_error(wi,"Too many parameters");
633 for(i=nn; (i<nrfp); i++)
637 for(i=0; (i<nral); i++) {
638 if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
639 gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
640 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
641 gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
643 for(i=0; (i<nrfp); i++)
645 push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line,wi);
649 void push_dihedraltype(directive d,t_params bt[],
650 t_bond_atomtype bat,char *line,
653 const char *formal[MAXATOMLIST+1] = {
662 const char *formnl[MAXATOMLIST+1] = {
668 "%*s%*s%*s%*s%*s%*s",
669 "%*s%*s%*s%*s%*s%*s%*s"
671 const char *formlf[MAXFORCEPARAM] = {
677 "%lf%lf%lf%lf%lf%lf",
678 "%lf%lf%lf%lf%lf%lf%lf",
679 "%lf%lf%lf%lf%lf%lf%lf%lf",
680 "%lf%lf%lf%lf%lf%lf%lf%lf%lf",
681 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
682 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
683 "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
685 int i,ft,ftype,nn,nrfp,nrfpA,nrfpB,nral;
687 char alc[MAXATOMLIST+1][20];
688 double c[MAXFORCEPARAM];
690 gmx_bool bAllowRepeat;
693 /* This routine accepts dihedraltypes defined from either 2 or 4 atoms.
695 * We first check for 2 atoms with the 3th column being an integer
696 * defining the type. If this isn't the case, we try it with 4 atoms
697 * and the 5th column defining the dihedral type.
699 nn=sscanf(line,formal[4],alc[0],alc[1],alc[2],alc[3],alc[4]);
700 if(nn>=3 && strlen(alc[2])==1 && isdigit(alc[2][0])) {
702 ft = strtol(alc[nral],NULL,10);
703 /* Move atom types around a bit and use 'X' for wildcard atoms
704 * to create a 4-atom dihedral definition with arbitrary atoms in
708 /* improper - the two atomtypes are 1,4. Use wildcards for 2,3 */
709 strcpy(alc[3],alc[1]);
712 /* alc[0] stays put */
714 /* proper - the two atomtypes are 2,3. Use wildcards for 1,4 */
716 strcpy(alc[2],alc[1]);
717 strcpy(alc[1],alc[0]);
720 } else if(nn==5 && strlen(alc[4])==1 && isdigit(alc[4][0])) {
722 ft = strtol(alc[nral],NULL,10);
724 sprintf(errbuf,"Incorrect number of atomtypes for dihedral (%d instead of 2 or 4)",nn-1);
725 warning_error(wi,errbuf);
731 /* Previously, we have always overwritten parameters if e.g. a torsion
732 with the same atomtypes occurs on multiple lines. However, CHARMM and
733 some other force fields specify multiple dihedrals over some bonds,
734 including cosines with multiplicity 6 and somethimes even higher.
735 Thus, they cannot be represented with Ryckaert-Bellemans terms.
736 To add support for these force fields, Dihedral type 9 is identical to
737 normal proper dihedrals, but repeated entries are allowed.
744 bAllowRepeat = FALSE;
748 ftype = ifunc_index(d,ft);
750 nrfpA = interaction_function[ftype].nrfpA;
751 nrfpB = interaction_function[ftype].nrfpB;
753 strcpy(f1,formnl[nral]);
754 strcat(f1,formlf[nrfp-1]);
756 /* Check number of parameters given */
757 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]))
760 /* Copy the B-state from the A-state */
761 copy_B_from_A(ftype,c);
764 warning_error(wi,"Not enough parameters");
765 } else if (nn > nrfpA && nn < nrfp) {
766 warning_error(wi,"Too many parameters or not enough parameters for topology B");
767 } else if (nn > nrfp) {
768 warning_error(wi,"Too many parameters");
770 for(i=nn; (i<nrfp); i++)
775 for(i=0; (i<4); i++) {
776 if(!strcmp(alc[i],"X"))
779 if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
780 gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
783 for(i=0; (i<nrfp); i++)
785 /* Always use 4 atoms here, since we created two wildcard atoms
786 * if there wasn't of them 4 already.
788 push_bondtype (&(bt[ftype]),&p,4,ftype,bAllowRepeat,line,wi);
792 void push_nbt(directive d,t_nbparam **nbt,gpp_atomtype_t atype,
793 char *pline,int nb_funct,
797 const char *form2="%*s%*s%*s%lf%lf";
798 const char *form3="%*s%*s%*s%lf%lf%lf";
799 const char *form4="%*s%*s%*s%lf%lf%lf%lf";
800 const char *form5="%*s%*s%*s%lf%lf%lf%lf%lf";
802 int i,f,n,ftype,atnr,nrfp;
810 if (sscanf (pline,"%s%s%d",a0,a1,&f) != 3) {
815 ftype=ifunc_index(d,f);
817 if (ftype != nb_funct) {
818 sprintf(errbuf,"Trying to add %s while the default nonbond type is %s",
819 interaction_function[ftype].longname,
820 interaction_function[nb_funct].longname);
821 warning_error(wi,errbuf);
825 /* Get the force parameters */
827 if (ftype == F_LJ14) {
828 n = sscanf(pline,form4,&c[0],&c[1],&c[2],&c[3]);
833 /* When the B topology parameters are not set,
834 * copy them from topology A
836 for(i=n; i<nrfp; i++)
839 else if (ftype == F_LJC14_Q) {
840 n = sscanf(pline,form5,&c[0],&c[1],&c[2],&c[3],&dum);
842 incorrect_n_param(wi);
846 else if (nrfp == 2) {
847 if (sscanf(pline,form3,&c[0],&c[1],&dum) != 2) {
848 incorrect_n_param(wi);
852 else if (nrfp == 3) {
853 if (sscanf(pline,form4,&c[0],&c[1],&c[2],&dum) != 3) {
854 incorrect_n_param(wi);
859 gmx_fatal(FARGS,"Number of force parameters for nonbonded interactions is %d"
860 " in file %s, line %d",nrfp,__FILE__,__LINE__);
862 for(i=0; (i<nrfp); i++)
865 /* Put the parameters in the matrix */
866 if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
867 gmx_fatal(FARGS,"Atomtype %s not found",a0);
868 if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
869 gmx_fatal(FARGS,"Atomtype %s not found",a1);
870 nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
874 for (i=0; i<nrfp; i++)
875 bId = bId && (nbp->c[i] == cr[i]);
877 sprintf(errbuf,"Overriding non-bonded parameters,");
879 fprintf(stderr," old:");
880 for (i=0; i<nrfp; i++)
881 fprintf(stderr," %g",nbp->c[i]);
882 fprintf(stderr," new\n%s\n",pline);
886 for (i=0; i<nrfp; i++)
891 push_gb_params (gpp_atomtype_t at, char *line,
896 double radius,vol,surftens,gb_radius,S_hct;
897 char atypename[STRLEN];
900 if ( (nfield = sscanf(line,"%s%lf%lf%lf%lf%lf",atypename,&radius,&vol,&surftens,&gb_radius,&S_hct)) != 6)
902 sprintf(errbuf,"Too few gb parameters for type %s\n",atypename);
906 /* Search for atomtype */
907 atype = get_atomtype_type(atypename,at);
911 printf("Couldn't find topology match for atomtype %s\n",atypename);
915 set_atomtype_gbparam(at,atype,radius,vol,surftens,gb_radius,S_hct);
919 push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
920 t_bond_atomtype bat, char *line,
923 const char *formal="%s%s%s%s%s%s%s%s";
925 int i,j,ft,ftype,nn,nrfp,nrfpA,nrfpB;
927 int nxcmap,nycmap,ncmap,read_cmap,sl,nct;
928 char s[20],alc[MAXATOMLIST+2][20];
930 gmx_bool bAllowRepeat;
933 /* Keep the compiler happy */
937 if((nn=sscanf(line,formal,alc[0],alc[1],alc[2],alc[3],alc[4],alc[5],alc[6],alc[7])) != nral+3)
939 sprintf(errbuf,"Incorrect number of atomtypes for cmap (%d instead of 5)",nn-1);
940 warning_error(wi,errbuf);
944 /* Compute an offset for each line where the cmap parameters start
945 * ie. where the atom types and grid spacing information ends
949 start += (int)strlen(alc[i]);
952 /* There are nn-1 spaces between the atom types and the grid spacing info in the cmap.itp file */
953 /* start is the position on the line where we start to read the actual cmap grid data from the itp file */
954 start = start + nn -1;
956 ft = strtol(alc[nral],NULL,10);
957 nxcmap = strtol(alc[nral+1],NULL,10);
958 nycmap = strtol(alc[nral+2],NULL,10);
960 /* Check for equal grid spacing in x and y dims */
963 gmx_fatal(FARGS,"Not the same grid spacing in x and y for cmap grid: x=%d, y=%d",nxcmap,nycmap);
966 ncmap = nxcmap*nycmap;
967 ftype = ifunc_index(d,ft);
968 nrfpA = strtol(alc[6],NULL,10)*strtol(alc[6],NULL,10);
969 nrfpB = strtol(alc[7],NULL,10)*strtol(alc[7],NULL,10);
972 /* Allocate memory for the CMAP grid */
974 srenew(bt->cmap,bt->ncmap);
976 /* Read in CMAP parameters */
980 while(isspace(*(line+start+sl)))
984 nn=sscanf(line+start+sl," %s ",s);
986 bt->cmap[i+(bt->ncmap)-nrfp]=strtod(s,NULL);
994 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]);
999 /* Check do that we got the number of parameters we expected */
1000 if(read_cmap==nrfpA)
1002 for(i=0;i<ncmap;i++)
1004 bt->cmap[i+ncmap]=bt->cmap[i];
1011 warning_error(wi,"Not enough cmap parameters");
1013 else if(read_cmap > nrfpA && read_cmap < nrfp)
1015 warning_error(wi,"Too many cmap parameters or not enough parameters for topology B");
1017 else if(read_cmap>nrfp)
1019 warning_error(wi,"Too many cmap parameters");
1024 /* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
1025 * so we can safely assign them each time
1027 bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
1028 bt->nc = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
1029 nct = (nral+1) * bt->nc;
1031 /* Allocate memory for the cmap_types information */
1032 srenew(bt->cmap_types,nct);
1034 for(i=0; (i<nral); i++)
1036 if(at && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1038 gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
1040 else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
1042 gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
1045 /* Assign a grid number to each cmap_type */
1046 bt->cmap_types[bt->nct++]=get_bond_atomtype_type(alc[i],bat);
1049 /* Assign a type number to this cmap */
1050 bt->cmap_types[bt->nct++]=bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
1052 /* Check for the correct number of atoms (again) */
1055 gmx_fatal(FARGS,"Incorrect number of atom types (%d) in cmap type %d\n",nct,bt->nc);
1058 /* Is this correct?? */
1059 for(i=0;i<MAXFORCEPARAM;i++)
1064 /* Push the bond to the bondlist */
1065 push_bondtype (&(bt[ftype]),&p,nral,ftype,FALSE,line,wi);
1069 static void push_atom_now(t_symtab *symtab,t_atoms *at,int atomnr,
1071 int type,char *ctype,int ptype,
1072 char *resnumberic,int cgnumber,
1073 char *resname,char *name,real m0,real q0,
1074 int typeB,char *ctypeB,real mB,real qB)
1076 int j,resind=0,resnr;
1080 if (((nr==0) && (atomnr != 1)) || (nr && (atomnr != at->nr+1)))
1081 gmx_fatal(FARGS,"Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = %d, while at->nr = %d)",atomnr,at->nr);
1083 j = strlen(resnumberic) - 1;
1084 if (isdigit(resnumberic[j])) {
1087 ric = resnumberic[j];
1088 if (j == 0 || !isdigit(resnumberic[j-1])) {
1089 gmx_fatal(FARGS,"Invalid residue number '%s' for atom %d",
1090 resnumberic,atomnr);
1093 resnr = strtol(resnumberic,NULL,10);
1096 resind = at->atom[nr-1].resind;
1098 if (nr == 0 || strcmp(resname,*at->resinfo[resind].name) != 0 ||
1099 resnr != at->resinfo[resind].nr ||
1100 ric != at->resinfo[resind].ic) {
1106 at->nres = resind + 1;
1107 srenew(at->resinfo,at->nres);
1108 at->resinfo[resind].name = put_symtab(symtab,resname);
1109 at->resinfo[resind].nr = resnr;
1110 at->resinfo[resind].ic = ric;
1112 resind = at->atom[at->nr-1].resind;
1115 /* New atom instance
1116 * get new space for arrays
1118 srenew(at->atom,nr+1);
1119 srenew(at->atomname,nr+1);
1120 srenew(at->atomtype,nr+1);
1121 srenew(at->atomtypeB,nr+1);
1124 at->atom[nr].type = type;
1125 at->atom[nr].ptype = ptype;
1126 at->atom[nr].q = q0;
1127 at->atom[nr].m = m0;
1128 at->atom[nr].typeB = typeB;
1129 at->atom[nr].qB = qB;
1130 at->atom[nr].mB = mB;
1132 at->atom[nr].resind = resind;
1133 at->atom[nr].atomnumber = atomicnumber;
1134 at->atomname[nr] = put_symtab(symtab,name);
1135 at->atomtype[nr] = put_symtab(symtab,ctype);
1136 at->atomtypeB[nr] = put_symtab(symtab,ctypeB);
1140 void push_cg(t_block *block, int *lastindex, int index, int a)
1143 fprintf (debug,"Index %d, Atom %d\n",index,a);
1145 if (((block->nr) && (*lastindex != index)) || (!block->nr)) {
1146 /* add a new block */
1148 srenew(block->index,block->nr+1);
1150 block->index[block->nr] = a + 1;
1154 void push_atom(t_symtab *symtab,t_block *cgs,
1155 t_atoms *at,gpp_atomtype_t atype,char *line,int *lastcg,
1159 int cgnumber,atomnr,type,typeB,nscan;
1160 char id[STRLEN],ctype[STRLEN],ctypeB[STRLEN],
1161 resnumberic[STRLEN],resname[STRLEN],name[STRLEN],check[STRLEN];
1165 /* Make a shortcut for writing in this molecule */
1168 /* Fixed parameters */
1169 if (sscanf(line,"%s%s%s%s%s%d",
1170 id,ctype,resnumberic,resname,name,&cgnumber) != 6) {
1174 sscanf(id,"%d",&atomnr);
1175 if ((type = get_atomtype_type(ctype,atype)) == NOTSET)
1176 gmx_fatal(FARGS,"Atomtype %s not found",ctype);
1177 ptype = get_atomtype_ptype(type,atype);
1179 /* Set default from type */
1180 q0 = get_atomtype_qA(type,atype);
1181 m0 = get_atomtype_massA(type,atype);
1186 /* Optional parameters */
1187 nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
1188 &q,&m,ctypeB,&qb,&mb,check);
1190 /* Nasty switch that falls thru all the way down! */
1196 if ((typeB = get_atomtype_type(ctypeB,atype)) == NOTSET) {
1197 gmx_fatal(FARGS,"Atomtype %s not found",ctypeB);
1199 qB = get_atomtype_qA(typeB,atype);
1200 mB = get_atomtype_massA(typeB,atype);
1206 warning_error(wi,"Too many parameters");
1213 fprintf(debug,"mB=%g, qB=%g, typeB=%d\n",mB,qB,typeB);
1215 push_cg(cgs,lastcg,cgnumber,nr);
1217 push_atom_now(symtab,at,atomnr,get_atomtype_atomnumber(type,atype),
1218 type,ctype,ptype,resnumberic,cgnumber,
1219 resname,name,m0,q0,typeB,
1220 typeB==type ? ctype : ctypeB,mB,qB);
1223 void push_molt(t_symtab *symtab,int *nmol,t_molinfo **mol,char *line,
1230 if ((sscanf(line,"%s%d",type,&nrexcl)) != 2) {
1231 warning_error(wi,"Expected a molecule type name and nrexcl");
1234 /* Test if this atomtype overwrites another */
1237 if (gmx_strcasecmp(*((*mol)[i].name),type) == 0)
1238 gmx_fatal(FARGS,"moleculetype %s is redefined",type);
1244 newmol = &((*mol)[*nmol-1]);
1245 init_molinfo(newmol);
1247 /* Fill in the values */
1248 newmol->name = put_symtab(symtab,type);
1249 newmol->nrexcl = nrexcl;
1250 newmol->excl_set = FALSE;
1253 static gmx_bool default_nb_params(int ftype,t_params bt[],t_atoms *at,
1254 t_param *p,int c_start,gmx_bool bB,gmx_bool bGenPairs)
1256 int i,j,ti,tj,ntype;
1259 int nr = bt[ftype].nr;
1260 int nral = NRAL(ftype);
1261 int nrfp = interaction_function[ftype].nrfpA;
1262 int nrfpB= interaction_function[ftype].nrfpB;
1264 if ((!bB && nrfp == 0) || (bB && nrfpB == 0))
1269 /* First test the generated-pair position to save
1270 * time when we have 1000*1000 entries for e.g. OPLS...
1274 ti=at->atom[p->a[0]].typeB;
1275 tj=at->atom[p->a[1]].typeB;
1277 ti=at->atom[p->a[0]].type;
1278 tj=at->atom[p->a[1]].type;
1280 pi=&(bt[ftype].param[ntype*ti+tj]);
1281 bFound=((ti == pi->a[0]) && (tj == pi->a[1]));
1284 /* Search explicitly if we didnt find it */
1286 for (i=0; ((i < nr) && !bFound); i++) {
1287 pi=&(bt[ftype].param[i]);
1289 for (j=0; ((j < nral) &&
1290 (at->atom[p->a[j]].typeB == pi->a[j])); j++);
1292 for (j=0; ((j < nral) &&
1293 (at->atom[p->a[j]].type == pi->a[j])); j++);
1300 if (nrfp+nrfpB > MAXFORCEPARAM) {
1301 gmx_incons("Too many force parameters");
1303 for (j=c_start; (j < nrfpB); j++)
1304 p->c[nrfp+j] = pi->c[j];
1307 for (j=c_start; (j < nrfp); j++)
1311 for (j=c_start; (j < nrfp); j++)
1317 static gmx_bool default_cmap_params(int ftype, t_params bondtype[],
1318 t_atoms *at, gpp_atomtype_t atype,
1319 t_param *p, gmx_bool bB,
1320 int *cmap_type, int *nparam_def)
1322 int i,j,nparam_found;
1324 gmx_bool bFound=FALSE;
1329 /* Match the current cmap angle against the list of cmap_types */
1330 for(i=0;i<bondtype->nct && !bFound;i+=6)
1339 (get_atomtype_batype(at->atom[p->a[0]].type,atype)==bondtype->cmap_types[i]) &&
1340 (get_atomtype_batype(at->atom[p->a[1]].type,atype)==bondtype->cmap_types[i+1]) &&
1341 (get_atomtype_batype(at->atom[p->a[2]].type,atype)==bondtype->cmap_types[i+2]) &&
1342 (get_atomtype_batype(at->atom[p->a[3]].type,atype)==bondtype->cmap_types[i+3]) &&
1343 (get_atomtype_batype(at->atom[p->a[4]].type,atype)==bondtype->cmap_types[i+4]))
1345 /* Found cmap torsion */
1347 ct = bondtype->cmap_types[i+5];
1353 /* If we did not find a matching type for this cmap torsion */
1356 gmx_fatal(FARGS,"Unknown cmap torsion between atoms %d %d %d %d %d\n",
1357 p->a[0]+1,p->a[1]+1,p->a[2]+1,p->a[3]+1,p->a[4]+1);
1360 *nparam_def = nparam_found;
1366 static gmx_bool default_params(int ftype,t_params bt[],
1367 t_atoms *at,gpp_atomtype_t atype,
1368 t_param *p,gmx_bool bB,
1369 t_param **param_def,
1372 int i,j,nparam_found;
1373 gmx_bool bFound,bSame;
1376 int nr = bt[ftype].nr;
1377 int nral = NRAL(ftype);
1378 int nrfpA= interaction_function[ftype].nrfpA;
1379 int nrfpB= interaction_function[ftype].nrfpB;
1381 if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
1385 /* We allow wildcards now. The first type (with or without wildcards) that
1386 * fits is used, so you should probably put the wildcarded bondtypes
1387 * at the end of each section.
1391 /* OPLS uses 1000s of dihedraltypes, so in order to speed up the scanning we have a
1392 * special case for this. Check for B state outside loop to speed it up.
1394 if( ftype==F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS)
1398 for (i=0; ((i < nr) && !bFound); i++)
1400 pi=&(bt[ftype].param[i]);
1403 ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
1404 ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
1405 ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
1406 ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
1413 for (i=0; ((i < nr) && !bFound); i++)
1415 pi=&(bt[ftype].param[i]);
1418 ((pi->AI==-1) || (get_atomtype_batype(at->atom[p->AI].type,atype)==pi->AI)) &&
1419 ((pi->AJ==-1) || (get_atomtype_batype(at->atom[p->AJ].type,atype)==pi->AJ)) &&
1420 ((pi->AK==-1) || (get_atomtype_batype(at->atom[p->AK].type,atype)==pi->AK)) &&
1421 ((pi->AL==-1) || (get_atomtype_batype(at->atom[p->AL].type,atype)==pi->AL))
1425 /* Find additional matches for this dihedral - necessary for ftype==9 which is used e.g. for charmm.
1426 * The rules in that case is that additional matches HAVE to be on adjacent lines!
1432 /* Continue from current i value */
1433 for(j=i+1 ; j<nr && bSame ; j+=2)
1435 pj=&(bt[ftype].param[j]);
1436 bSame = (pi->AI == pj->AI && pi->AJ == pj->AJ && pi->AK == pj->AK && pi->AL == pj->AL);
1439 /* nparam_found will be increased as long as the numbers match */
1442 } else { /* Not a dihedral */
1443 for (i=0; ((i < nr) && !bFound); i++) {
1444 pi=&(bt[ftype].param[i]);
1446 for (j=0; ((j < nral) &&
1447 (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
1450 for (j=0; ((j < nral) &&
1451 (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
1460 *nparam_def = nparam_found;
1467 void push_bond(directive d,t_params bondtype[],t_params bond[],
1468 t_atoms *at,gpp_atomtype_t atype,char *line,
1469 gmx_bool bBonded,gmx_bool bGenPairs,real fudgeQQ,
1470 gmx_bool bZero,gmx_bool *bWarn_copy_A_B,
1473 const char *aaformat[MAXATOMLIST]= {
1481 const char *asformat[MAXATOMLIST]= {
1486 "%*s%*s%*s%*s%*s%*s",
1487 "%*s%*s%*s%*s%*s%*s%*s"
1489 const char *ccformat="%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf";
1490 int nr,i,j,nral,nread,ftype;
1491 char format[STRLEN];
1492 /* One force parameter more, so we can check if we read too many */
1493 double cc[MAXFORCEPARAM+1];
1494 int aa[MAXATOMLIST+1];
1495 t_param param,paramB,*param_defA,*param_defB;
1496 gmx_bool bFoundA=FALSE,bFoundB=FALSE,bDef,bPert,bSwapParity=FALSE;
1497 int nparam_defA,nparam_defB;
1500 nparam_defA=nparam_defB=0;
1502 ftype = ifunc_index(d,1);
1504 for(j=0; j<MAXATOMLIST; j++)
1506 bDef = (NRFP(ftype) > 0);
1508 nread = sscanf(line,aaformat[nral-1],
1509 &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1513 } else if (nread == nral)
1514 ftype = ifunc_index(d,1);
1516 /* this is a hack to allow for virtual sites with swapped parity */
1517 bSwapParity = (aa[nral]<0);
1519 aa[nral] = -aa[nral];
1520 ftype = ifunc_index(d,aa[nral]);
1527 gmx_fatal(FARGS,"Negative function types only allowed for %s and %s",
1528 interaction_function[F_VSITE3FAD].longname,
1529 interaction_function[F_VSITE3OUT].longname);
1533 /* Check for double atoms and atoms out of bounds */
1534 for(i=0; (i<nral); i++) {
1535 if ( aa[i] < 1 || aa[i] > at->nr )
1536 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1537 "Atom index (%d) in %s out of bounds (1-%d).\n"
1538 "This probably means that you have inserted topology section \"%s\"\n"
1539 "in a part belonging to a different molecule than you intended to.\n"
1540 "In that case move the \"%s\" section to the right molecule.",
1541 get_warning_file(wi),get_warning_line(wi),
1542 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1543 for(j=i+1; (j<nral); j++)
1544 if (aa[i] == aa[j]) {
1545 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1549 if (ftype == F_SETTLE)
1550 if (aa[0]+2 > at->nr)
1551 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1552 " Atom index (%d) in %s out of bounds (1-%d)\n"
1553 " Settle works on atoms %d, %d and %d",
1554 get_warning_file(wi),get_warning_line(wi),
1555 aa[0],dir2str(d),at->nr,
1556 aa[0],aa[0]+1,aa[0]+2);
1558 /* default force parameters */
1559 for(j=0; (j<MAXATOMLIST); j++)
1560 param.a[j] = aa[j]-1;
1561 for(j=0; (j<MAXFORCEPARAM); j++)
1564 /* Get force params for normal and free energy perturbation
1565 * studies, as determined by types!
1569 bFoundA = default_params(ftype,bondtype,at,atype,¶m,FALSE,¶m_defA,&nparam_defA);
1571 /* Copy the A-state and B-state default parameters */
1572 for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1573 param.c[j] = param_defA->c[j];
1575 bFoundB = default_params(ftype,bondtype,at,atype,¶m,TRUE,¶m_defB,&nparam_defB);
1577 /* Copy only the B-state default parameters */
1578 for(j=NRFPA(ftype); (j<NRFP(ftype)); j++)
1579 param.c[j] = param_defB->c[j];
1581 } else if (ftype == F_LJ14) {
1582 bFoundA = default_nb_params(ftype, bondtype,at,¶m,0,FALSE,bGenPairs);
1583 bFoundB = default_nb_params(ftype, bondtype,at,¶m,0,TRUE, bGenPairs);
1584 } else if (ftype == F_LJC14_Q) {
1585 param.c[0] = fudgeQQ;
1586 /* Fill in the A-state charges as default parameters */
1587 param.c[1] = at->atom[param.a[0]].q;
1588 param.c[2] = at->atom[param.a[1]].q;
1589 /* The default LJ parameters are the standard 1-4 parameters */
1590 bFoundA = default_nb_params(F_LJ14,bondtype,at,¶m,3,FALSE,bGenPairs);
1592 } else if (ftype == F_LJC_PAIRS_NB) {
1593 /* Defaults are not supported here */
1597 gmx_incons("Unknown function type in push_bond");
1601 /* Manually specified parameters - in this case we discard multiple torsion info! */
1603 strcpy(format,asformat[nral-1]);
1604 strcat(format,ccformat);
1606 nread = sscanf(line,format,&cc[0],&cc[1],&cc[2],&cc[3],&cc[4],&cc[5],
1607 &cc[6],&cc[7],&cc[8],&cc[9],&cc[10],&cc[11],&cc[12]);
1609 if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0))
1611 /* We only have to issue a warning if these atoms are perturbed! */
1613 for(j=0; (j<nral); j++)
1614 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1616 if (bPert && *bWarn_copy_A_B)
1619 "Some parameters for bonded interaction involving perturbed atoms are specified explicitly in state A, but not B - copying A to B");
1621 *bWarn_copy_A_B = FALSE;
1624 /* If only the A parameters were specified, copy them to the B state */
1625 /* The B-state parameters correspond to the first nrfpB
1626 * A-state parameters.
1628 for(j=0; (j<NRFPB(ftype)); j++)
1629 cc[nread++] = cc[j];
1632 /* If nread was 0 or EOF, no parameters were read => use defaults.
1633 * If nread was nrfpA we copied above so nread=nrfp.
1634 * If nread was nrfp we are cool.
1635 * For F_LJC14_Q we allow supplying fudgeQQ only.
1636 * Anything else is an error!
1638 if ((nread != 0) && (nread != EOF) && (nread != NRFP(ftype)) &&
1639 !(ftype == F_LJC14_Q && nread == 1))
1641 gmx_fatal(FARGS,"Incorrect number of parameters - found %d, expected %d or %d for %s.",
1642 nread,NRFPA(ftype),NRFP(ftype),
1643 interaction_function[ftype].longname);
1646 for(j=0; (j<nread); j++)
1649 /* Check whether we have to use the defaults */
1650 if (nread == NRFP(ftype))
1657 /* nread now holds the number of force parameters read! */
1662 /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */
1665 if ((nparam_defA != nparam_defB) || ((nparam_defA>1 || nparam_defB>1) && (param_defA!=param_defB)))
1668 "Cannot automatically perturb a torsion with multiple terms to different form.\n"
1669 "Please specify perturbed parameters manually for this torsion in your topology!");
1670 warning_error(wi,errbuf);
1674 if (nread > 0 && nread < NRFPA(ftype))
1676 /* Issue an error, do not use defaults */
1677 sprintf(errbuf,"Not enough parameters, there should be at least %d (or 0 for defaults)",NRFPA(ftype));
1678 warning_error(wi,errbuf);
1681 if (nread == 0 || nread == EOF)
1685 if (interaction_function[ftype].flags & IF_VSITE)
1687 /* set them to NOTSET, will be calculated later */
1688 for(j=0; (j<MAXFORCEPARAM); j++)
1689 param.c[j] = NOTSET;
1692 param.C1 = -1; /* flag to swap parity of vsite construction */
1698 fprintf(stderr,"NOTE: No default %s types, using zeroes\n",
1699 interaction_function[ftype].longname);
1703 sprintf(errbuf,"No default %s types",interaction_function[ftype].longname);
1704 warning_error(wi,errbuf);
1715 param.C0 = 360 - param.C0;
1718 param.C2 = -param.C2;
1725 /* We only have to issue a warning if these atoms are perturbed! */
1727 for(j=0; (j<nral); j++)
1728 bPert = bPert || PERTURBED(at->atom[param.a[j]]);
1732 sprintf(errbuf,"No default %s types for perturbed atoms, "
1733 "using normal values",interaction_function[ftype].longname);
1740 if ((ftype==F_PDIHS || ftype==F_ANGRES || ftype==F_ANGRESZ)
1741 && param.c[5]!=param.c[2])
1742 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1743 " %s multiplicity can not be perturbed %f!=%f",
1744 get_warning_file(wi),get_warning_line(wi),
1745 interaction_function[ftype].longname,
1746 param.c[2],param.c[5]);
1748 if (IS_TABULATED(ftype) && param.c[2]!=param.c[0])
1749 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1750 " %s table number can not be perturbed %d!=%d",
1751 get_warning_file(wi),get_warning_line(wi),
1752 interaction_function[ftype].longname,
1753 (int)(param.c[0]+0.5),(int)(param.c[2]+0.5));
1755 /* Dont add R-B dihedrals where all parameters are zero (no interaction) */
1756 if (ftype==F_RBDIHS) {
1758 for(i=0;i<NRFP(ftype);i++) {
1766 /* Put the values in the appropriate arrays */
1767 add_param_to_list (&bond[ftype],¶m);
1769 /* Push additional torsions from FF for ftype==9 if we have them.
1770 * We have already checked that the A/B states do not differ in this case,
1771 * so we do not have to double-check that again, or the vsite stuff.
1772 * In addition, those torsions cannot be automatically perturbed.
1774 if(bDef && ftype==F_PDIHS)
1776 for(i=1;i<nparam_defA;i++)
1778 /* Advance pointer! */
1780 for(j=0; (j<NRFPA(ftype)+NRFPB(ftype)); j++)
1781 param.c[j] = param_defA->c[j];
1782 /* And push the next term for this torsion */
1783 add_param_to_list (&bond[ftype],¶m);
1788 void push_cmap(directive d, t_params bondtype[], t_params bond[],
1789 t_atoms *at, gpp_atomtype_t atype, char *line,
1790 gmx_bool *bWarn_copy_A_B,
1793 const char *aaformat[MAXATOMLIST+1]=
1804 int i,j,nr,ftype,nral,nread,ncmap_params;
1806 int aa[MAXATOMLIST+1];
1809 t_param param,paramB,*param_defA,*param_defB;
1811 ftype = ifunc_index(d,1);
1813 nr = bondtype[ftype].nr;
1816 nread = sscanf(line,aaformat[nral-1],
1817 &aa[0],&aa[1],&aa[2],&aa[3],&aa[4],&aa[5]);
1824 else if (nread == nral)
1826 ftype = ifunc_index(d,1);
1829 /* Check for double atoms and atoms out of bounds */
1832 if(aa[i] < 1 || aa[i] > at->nr)
1834 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1835 "Atom index (%d) in %s out of bounds (1-%d).\n"
1836 "This probably means that you have inserted topology section \"%s\"\n"
1837 "in a part belonging to a different molecule than you intended to.\n"
1838 "In that case move the \"%s\" section to the right molecule.",
1839 get_warning_file(wi),get_warning_line(wi),
1840 aa[i],dir2str(d),at->nr,dir2str(d),dir2str(d));
1843 for(j=i+1; (j<nral); j++)
1847 sprintf(errbuf,"Duplicate atom index (%d) in %s",aa[i],dir2str(d));
1853 /* default force parameters */
1854 for(j=0; (j<MAXATOMLIST); j++)
1856 param.a[j] = aa[j]-1;
1858 for(j=0; (j<MAXFORCEPARAM); j++)
1863 /* Get the cmap type for this cmap angle */
1864 bFound = default_cmap_params(ftype,bondtype,at,atype,¶m,FALSE,&cmap_type,&ncmap_params);
1866 /* We want exactly one parameter (the cmap type in state A (currently no state B) back */
1867 if(bFound && ncmap_params==1)
1869 /* Put the values in the appropriate arrays */
1870 param.c[0]=cmap_type;
1871 add_param_to_list(&bond[ftype],¶m);
1875 /* This is essentially the same check as in default_cmap_params() done one more time */
1876 gmx_fatal(FARGS, "Unable to assign a cmap type to torsion %d %d %d %d and %d\n",
1877 param.a[0]+1,param.a[1]+1,param.a[2]+1,param.a[3]+1,param.a[4]+1);
1883 void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
1884 t_atoms *at,gpp_atomtype_t atype,char *line,
1888 int type,ftype,j,n,ret,nj,a;
1890 double *weight=NULL,weight_tot;
1893 /* default force parameters */
1894 for(j=0; (j<MAXATOMLIST); j++)
1895 param.a[j] = NOTSET;
1896 for(j=0; (j<MAXFORCEPARAM); j++)
1900 ret = sscanf(ptr,"%d%n",&a,&n);
1903 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1904 " Expected an atom index in section \"%s\"",
1905 get_warning_file(wi),get_warning_line(wi),
1910 ret = sscanf(ptr,"%d%n",&type,&n);
1912 ftype = ifunc_index(d,type);
1917 ret = sscanf(ptr,"%d%n",&a,&n);
1922 srenew(weight,nj+20);
1930 /* Here we use the A-state mass as a parameter.
1931 * Note that the B-state mass has no influence.
1933 weight[nj] = at->atom[atc[nj]].m;
1937 ret = sscanf(ptr,"%lf%n",&(weight[nj]),&n);
1940 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1941 " No weight or negative weight found for vsiten constructing atom %d (atom index %d)",
1942 get_warning_file(wi),get_warning_line(wi),
1946 gmx_fatal(FARGS,"Unknown vsiten type %d",type);
1948 weight_tot += weight[nj];
1954 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1955 " Expected more than one atom index in section \"%s\"",
1956 get_warning_file(wi),get_warning_line(wi),
1959 if (weight_tot == 0)
1960 gmx_fatal(FARGS,"[ file %s, line %d ]:\n"
1961 " The total mass of the construting atoms is zero",
1962 get_warning_file(wi),get_warning_line(wi));
1964 for(j=0; j<nj; j++) {
1965 param.a[1] = atc[j];
1967 param.c[1] = weight[j]/weight_tot;
1968 /* Put the values in the appropriate arrays */
1969 add_param_to_list (&bond[ftype],¶m);
1976 void push_mol(int nrmols,t_molinfo mols[],char *pline,int *whichmol,
1984 if (sscanf(pline,"%s%d",type,&copies) != 2) {
1989 /* search moleculename */
1990 for (i=0; ((i<nrmols) && gmx_strcasecmp(type,*(mols[i].name))); i++)
1997 gmx_fatal(FARGS,"No such moleculetype %s",type);
2000 void init_block2(t_block2 *b2, int natom)
2005 snew(b2->nra,b2->nr);
2007 for(i=0; (i<b2->nr); i++)
2011 void done_block2(t_block2 *b2)
2016 for(i=0; (i<b2->nr); i++)
2024 void push_excl(char *line, t_block2 *b2)
2028 char base[STRLEN],format[STRLEN];
2030 if (sscanf(line,"%d",&i)==0)
2033 if ((1 <= i) && (i <= b2->nr))
2037 fprintf(debug,"Unbound atom %d\n",i-1);
2042 strcpy(format,base);
2043 strcat(format,"%d");
2044 n=sscanf(line,format,&j);
2046 if ((1 <= j) && (j <= b2->nr)) {
2048 srenew(b2->a[i],++(b2->nra[i]));
2049 b2->a[i][b2->nra[i]-1]=j;
2050 /* also add the reverse exclusion! */
2051 srenew(b2->a[j],++(b2->nra[j]));
2052 b2->a[j][b2->nra[j]-1]=i;
2056 gmx_fatal(FARGS,"Invalid Atomnr j: %d, b2->nr: %d\n",j,b2->nr);
2061 void b_to_b2(t_blocka *b, t_block2 *b2)
2066 for(i=0; (i<b->nr); i++)
2067 for(j=b->index[i]; (j<b->index[i+1]); j++) {
2069 srenew(b2->a[i],++b2->nra[i]);
2070 b2->a[i][b2->nra[i]-1]=a;
2074 void b2_to_b(t_block2 *b2, t_blocka *b)
2080 for(i=0; (i<b2->nr); i++) {
2082 for(j=0; (j<b2->nra[i]); j++)
2083 b->a[nra+j]=b2->a[i][j];
2086 /* terminate list */
2090 static int icomp(const void *v1, const void *v2)
2092 return (*((atom_id *) v1))-(*((atom_id *) v2));
2095 void merge_excl(t_blocka *excl, t_block2 *b2)
2103 else if (b2->nr != excl->nr) {
2104 gmx_fatal(FARGS,"DEATH HORROR: b2->nr = %d, while excl->nr = %d",
2108 fprintf(debug,"Entering merge_excl\n");
2110 /* First copy all entries from excl to b2 */
2113 /* Count and sort the exclusions */
2115 for(i=0; (i<b2->nr); i++) {
2116 if (b2->nra[i] > 0) {
2117 /* remove double entries */
2118 qsort(b2->a[i],(size_t)b2->nra[i],(size_t)sizeof(b2->a[i][0]),icomp);
2120 for(j=1; (j<b2->nra[i]); j++)
2121 if (b2->a[i][j]!=b2->a[i][k-1]) {
2122 b2->a[i][k]=b2->a[i][j];
2130 srenew(excl->a,excl->nra);
2135 int add_atomtype_decoupled(t_symtab *symtab,gpp_atomtype_t at,
2136 t_nbparam ***nbparam,t_nbparam ***pair)
2142 /* Define an atom type with all parameters set to zero (no interactions) */
2145 /* Type for decoupled atoms could be anything,
2146 * this should be changed automatically later when required.
2148 atom.ptype = eptAtom;
2149 for (i=0; (i<MAXFORCEPARAM); i++)
2152 nr = add_atomtype(at,symtab,&atom,"decoupled",¶m,-1,0.0,0.0,0.0,0,0,0);
2154 /* Add space in the non-bonded parameters matrix */
2155 realloc_nb_params(at,nbparam,pair);
2160 static void convert_pairs_to_pairsQ(t_params *plist,
2161 real fudgeQQ,t_atoms *atoms)
2167 /* Copy the pair list to the pairQ list */
2168 plist[F_LJC14_Q] = plist[F_LJ14];
2169 /* Empty the pair list */
2170 plist[F_LJ14].nr = 0;
2171 plist[F_LJ14].param = NULL;
2172 param = plist[F_LJC14_Q].param;
2173 for(i=0; i<plist[F_LJC14_Q].nr; i++) {
2176 param[i].c[0] = fudgeQQ;
2177 param[i].c[1] = atoms->atom[param[i].a[0]].q;
2178 param[i].c[2] = atoms->atom[param[i].a[1]].q;
2184 static void generate_LJCpairsNB(t_molinfo *mol,int nb_funct,t_params *nbp)
2193 atom = mol->atoms.atom;
2195 ntype = sqrt(nbp->nr);
2197 for (i=0; i<MAXATOMLIST; i++)
2198 param.a[i] = NOTSET;
2199 for (i=0; i<MAXFORCEPARAM; i++)
2200 param.c[i] = NOTSET;
2202 /* Add a pair interaction for all non-excluded atom pairs */
2204 for(i=0; i<n; i++) {
2205 for(j=i+1; j<n; j++) {
2207 for(k=excl->index[i]; k<excl->index[i+1]; k++) {
2208 if (excl->a[k] == j)
2212 if (nb_funct != F_LJ)
2213 gmx_fatal(FARGS,"Can only generate non-bonded pair interactions for Van der Waals type Lennard-Jones");
2216 param.c[0] = atom[i].q;
2217 param.c[1] = atom[j].q;
2218 param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0];
2219 param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1];
2220 add_param_to_list(&mol->plist[F_LJC_PAIRS_NB],¶m);
2226 static void set_excl_all(t_blocka *excl)
2230 /* Get rid of the current exclusions and exclude all atom pairs */
2232 excl->nra = nat*nat;
2233 srenew(excl->a,excl->nra);
2235 for(i=0; i<nat; i++) {
2237 for(j=0; j<nat; j++)
2240 excl->index[nat] = k;
2243 static void decouple_atoms(t_atoms *atoms,int atomtype_decouple,
2244 int couple_lam0,int couple_lam1)
2248 for(i=0; i<atoms->nr; i++) {
2249 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamVDW) {
2250 atoms->atom[i].q = 0.0;
2252 if (couple_lam0 == ecouplamNONE || couple_lam0 == ecouplamQ) {
2253 atoms->atom[i].type = atomtype_decouple;
2255 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamVDW) {
2256 atoms->atom[i].qB = 0.0;
2258 if (couple_lam1 == ecouplamNONE || couple_lam1 == ecouplamQ) {
2259 atoms->atom[i].typeB = atomtype_decouple;
2264 void convert_moltype_couple(t_molinfo *mol,int atomtype_decouple,real fudgeQQ,
2265 int couple_lam0,int couple_lam1,
2266 gmx_bool bCoupleIntra,int nb_funct,t_params *nbp)
2268 convert_pairs_to_pairsQ(mol->plist,fudgeQQ,&mol->atoms);
2269 if (!bCoupleIntra) {
2270 generate_LJCpairsNB(mol,nb_funct,nbp);
2271 set_excl_all(&mol->excls);
2273 decouple_atoms(&mol->atoms,atomtype_decouple,couple_lam0,couple_lam1);