t_block mols; /* Molecules */
t_blocka excls; /* Exclusions */
t_params plist[F_NRE]; /* Parameters in old style */
-
} t_molinfo;
-typedef struct {
- int nr; /* The number of atomtypes */
- t_atom *atom; /* Array of atoms */
- char ***atomname; /* Names of the atomtypes */
- t_param *nb; /* Nonbonded force default params */
- int *bondatomtype; /* The bond_atomtype for each atomtype */
- real *radius; /* Radius for GBSA stuff */
- real *vol; /* Effective volume for GBSA */
- real *surftens; /* Surface tension with water, for GBSA */
- int *atomnumber; /* Atomic number, used for QM/MM */
-} t_atomtype;
-
-typedef struct {
- int nr; /* Number of bond_atomtypes */
- char ***atomname; /* Names of the bond_atomtypes */
-} t_bond_atomtype;
-
-
typedef enum {
d_defaults,
d_atomtypes,
"invalid"
};
-extern void convert_harmonics(int nrmols,t_molinfo mols[],t_atomtype *atype);
-
#endif /* _grompp_h */
#include "grompp.h"
enum { eqgNone, eqgLinear, eqgYang, eqgBultinck,
- eqgSM1, eqgSM2, eqgSM3, eqgSM4, eqgNR };
+ eqgSMp, eqgSMs, eqgSMps, eqgSMg, eqgSMpg, eqgNR };
extern real generate_charges_sm(FILE *fp,char *molname,void *eem,
t_atoms *atoms,rvec x[],
rms += sqr(chi0-md->Chi0_0);
if (chi0 > md->Chi0_1)
rms += sqr(chi0-md->Chi0_1);
- if (md->eemtype != eqgSM1) {
+ if (md->eemtype != eqgSMp) {
w = gsl_vector_get(v, k++);
if (w < md->w_0)
rms += sqr(w-md->w_0);
my_func.f = &dipole_function;
my_func.n = md->nparam*2;
- if (md->eemtype != eqgSM1)
+ if (md->eemtype != eqgSMp)
my_func.n += md->nparam;
my_func.params = (void *) md;
gsl_vector_set (dx, k++, stepsize*chi0);
w = eem_get_w(md->eem,md->index[i]);
- if (md->eemtype != eqgSM1) {
+ if (md->eemtype != eqgSMp) {
w = guess_new_param(w,stepsize,md->w_0,md->w_1,rng,bRand);
gsl_vector_set (x, k, w);
gsl_vector_set (dx, k++, stepsize*w);
parse_common_args(&argc,argv,PCA_CAN_VIEW,NFILE,fnm,asize(pa),pa,
asize(desc),desc,0,NULL);
- eemtype = eqgSM1;
+ eemtype = eqgSMp;
if (qgen[0]) {
eemtype = name2eemtype(qgen[0]);
if (eemtype == -1)
- eemtype = eqgSM1;
+ eemtype = eqgSMp;
}
- if (eemtype > eqgSM3)
- gmx_fatal(FARGS,"Only models SM1 and SM2 implemented so far");
+ if (eemtype > eqgSMps)
+ gmx_fatal(FARGS,"Only models SMp, SMs and SMps implemented so far");
logf = fopen(opt2fn("-g",NFILE,fnm),"w");
md = read_moldip(logf,opt2fn("-f",NFILE,fnm),
gen_vsite.c gen_vsite.h \
genhydro.c genhydro.h \
gmxcpp.c gmxcpp.h \
+gpp_atomtype.c gpp_atomtype.h \
+gpp_bond_atomtype.c gpp_bond_atomtype.h \
h_db.c h_db.h \
hackblock.c hackblock.h \
hizzie.c hizzie.h \
sorting.c sorting.h \
specbond.c specbond.h \
ter_db.c ter_db.h \
-tomorse.c \
+tomorse.c gpp_tomorse.h \
topcat.c topcat.h \
topdirs.c topdirs.h \
topexcl.c topexcl.h \
#include "toputil.h"
#include "convparm.h"
#include "names.h"
+#include "gpp_atomtype.h"
static int round_check(real r,int limit,int ftype,char *name)
{
#include "index.h"
#include "names.h"
#include "futil.h"
+#include "gpp_atomtype.h"
#define MAXNAME 32
#define OPENDIR '[' /* starting sign for directive */
return type;
}
-static int nm2type(char *name, t_atomtype *atype)
+static int vsite_nm2type(char *name, t_atomtype atype)
{
- int tp,j;
+ int tp;
- tp=NOTSET;
- for(j=0; (j < atype->nr) && (tp==NOTSET); j++)
- if (strcasecmp(name,*(atype->atomname[j])) == 0)
- tp=j;
- if (tp==NOTSET)
- gmx_fatal(FARGS,"Dummy mass type (%s) not found in atom type database",name);
+ tp = get_atomtype_type(name,atype);
+ if (tp == NOTSET)
+ gmx_fatal(FARGS,"Dummy mass type (%s) not found in atom type database",
+ name);
+
return tp;
}
}
-static int gen_vsites_trp(t_atomtype *atype, rvec *newx[],
+static int gen_vsites_trp(t_atomtype atype, rvec *newx[],
t_atom *newatom[], char ***newatomname[],
int *o2n[], int *newvsite_type[], int *newcgnr[],
t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
}
/* get dummy mass type */
- tpM=nm2type("MW",atype);
+ tpM=vsite_nm2type("MW",atype);
/* make space for 2 masses: shift all atoms starting with CB */
i0=ats[atCB];
for(j=0; j<NMASS; j++)
}
-static int gen_vsites_tyr(t_atomtype *atype, rvec *newx[],
+static int gen_vsites_tyr(t_atomtype atype, rvec *newx[],
t_atom *newatom[], char ***newatomname[],
int *o2n[], int *newvsite_type[], int *newcgnr[],
t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
at->atom[ats[atHH]].m=at->atom[ats[atHH]].mB=0;
/* get dummy mass type */
- tpM=nm2type("MW",atype);
+ tpM=vsite_nm2type("MW",atype);
/* make space for 1 mass: shift HH only */
i0=ats[atHH];
atM=i0+*nadd;
static char atomnamesuffix[] = "1234";
-void do_vsites(int nrtp, t_restp rtp[], t_atomtype *atype,
+void do_vsites(int nrtp, t_restp rtp[], t_atomtype atype,
t_atoms *at, t_symtab *symtab, rvec *x[],
t_params plist[], int *vsite_type[], int *cgnr[],
real mHmult, bool bVsiteAromatics, char *ff)
&nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
/* get Heavy atom type */
tpHeavy=get_atype(Heavy,at,nrtp,rtp,aan);
- strcpy(tpname,type2nm(tpHeavy,atype));
+ strcpy(tpname,get_atomtype_name(tpHeavy,atype));
bWARNING=FALSE;
bAddVsiteParam=TRUE;
(*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
/* get dummy mass type from first char of heavy atom type (N or C) */
- strcpy(nexttpname,type2nm(get_atype(heavies[0],at,nrtp,rtp,aan),atype));
+ strcpy(nexttpname,get_atomtype_name(get_atype(heavies[0],at,nrtp,rtp,aan),atype));
ch=get_dummymass_name(vsiteconflist,nvsiteconf,tpname,nexttpname);
if(ch==NULL)
else
strcpy(name,ch);
- tpM=nm2type(name,atype);
+ tpM=vsite_nm2type(name,atype);
/* make space for 2 masses: shift all atoms starting with 'Heavy' */
#define NMASS 2
i0=Heavy;
#include "typedefs.h"
#include "grompp.h"
+#include "gpp_atomtype.h"
#include "hackblock.h"
/* stuff for pdb2gmx */
-extern void do_vsites(int nrtp, t_restp rtp[], t_atomtype *atype,
+extern void do_vsites(int nrtp, t_restp rtp[], t_atomtype atype,
t_atoms *at, t_symtab *symtab, rvec *x[],
t_params plist[], int *dummy_type[], int *cgnr[],
real mHmult, bool bVSiteAromatics, char *ff);
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include "smalloc.h"
+#include "sysstuff.h"
+#include "macros.h"
+#include "string2.h"
+#include "topdirs.h"
+#include "toputil.h"
+#include "topdirs.h"
+#include "toputil.h"
+#include "symtab.h"
+#include "gmx_fatal.h"
+#include "txtdump.h"
+#include "gpp_atomtype.h"
+
+typedef struct {
+ int nr; /* The number of atomtypes */
+ t_atom *atom; /* Array of atoms */
+ char ***atomname; /* Names of the atomtypes */
+ t_param *nb; /* Nonbonded force default params */
+ int *bondatomtype; /* The bond_atomtype for each atomtype */
+ real *radius; /* Radius for GBSA stuff */
+ real *vol; /* Effective volume for GBSA */
+ real *surftens; /* Surface tension with water, for GBSA */
+ int *atomnumber; /* Atomic number, used for QM/MM */
+} gpp_atomtype;
+
+int get_atomtype_type(char *str,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ int i;
+
+ for (i=0; (i<ga->nr); i++)
+ if (strcasecmp(str,*(ga->atomname[i])) == 0)
+ return i;
+
+ return NOTSET;
+}
+
+int get_atomtype_ntypes(t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ return ga->nr;
+}
+
+char *get_atomtype_name(int nt, t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NULL;
+
+ return *(ga->atomname[nt]);
+}
+
+real get_atomtype_massA(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atom[nt].m;
+}
+
+real get_atomtype_massB(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atom[nt].mB;
+}
+
+real get_atomtype_qA(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atom[nt].q;
+}
+
+real get_atomtype_qB(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atom[nt].qB;
+}
+
+int get_atomtype_ptype(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atom[nt].ptype;
+}
+
+int get_atomtype_batype(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->bondatomtype[nt];
+}
+
+int get_atomtype_atomnumber(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->atomnumber[nt];
+}
+
+real get_atomtype_radius(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->radius[nt];
+}
+
+real get_atomtype_vol(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->vol[nt];
+}
+
+real get_atomtype_surftens(int nt,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ return ga->surftens[nt];
+}
+
+real get_atomtype_nbparam(int nt,int param,t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+ if ((param < 0) || (param >= MAXFORCEPARAM))
+ return NOTSET;
+ return ga->nb[nt].c[param];
+}
+
+t_atomtype init_atomtype(void)
+{
+ gpp_atomtype *ga;
+
+ snew(ga,1);
+
+ return (t_atomtype ) ga;
+}
+
+int set_atomtype(int nt,t_atomtype at,t_symtab *tab,
+ t_atom *a,char *name,t_param *nb,
+ int bondatomtype,
+ real radius,real vol,real surftens,int atomnumber)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NOTSET;
+
+ ga->atom[nt] = *a;
+ ga->atomname[nt] = put_symtab(tab,name);
+ ga->nb[nt] = *nb;
+ ga->bondatomtype[nt] = bondatomtype;
+ ga->radius[nt] = radius;
+ ga->vol[nt] = vol;
+ ga->surftens[nt] = surftens;
+ ga->atomnumber[nt] = atomnumber;
+
+ return nt;
+}
+
+int add_atomtype(t_atomtype at,t_symtab *tab,
+ t_atom *a,char *name,t_param *nb,
+ int bondatomtype,
+ real radius,real vol,real surftens,real atomnumber)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ ga->nr++;
+ srenew(ga->atom,ga->nr);
+ srenew(ga->atomname,ga->nr);
+ srenew(ga->nb,ga->nr);
+ srenew(ga->bondatomtype,ga->nr);
+ srenew(ga->radius,ga->nr);
+ srenew(ga->vol,ga->nr);
+ srenew(ga->surftens,ga->nr);
+ srenew(ga->atomnumber,ga->nr);
+
+ return set_atomtype(ga->nr-1,at,tab,a,name,nb,bondatomtype,radius,
+ vol,surftens,atomnumber);
+}
+
+void print_at (FILE * out, t_atomtype at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ int i;
+ t_atom *atom = ga->atom;
+ t_param *nb = ga->nb;
+
+ fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
+ fprintf (out,"; %6s %8s %8s %8s %12s %12s\n",
+ "type","mass","charge","particle","c6","c12");
+ for (i=0; (i<ga->nr); i++)
+ fprintf(out,"%8s %8.3f %8.3f %8s %12e %12e\n",
+ *(ga->atomname[i]),atom[i].m,atom[i].q,"A",
+ nb[i].C0,nb[i].C1);
+
+ fprintf (out,"\n");
+}
+
+void done_atomtype(t_atomtype *at)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) *at;
+
+ sfree(ga->atom);
+ sfree(ga->atomname);
+ sfree(ga->nb);
+ sfree(ga->bondatomtype);
+ sfree(ga->radius);
+ sfree(ga->vol);
+ sfree(ga->surftens);
+ sfree(ga->atomnumber);
+ ga->nr = 0;
+ sfree(ga);
+ *at = NULL;
+}
+
+static int search_atomtypes(t_atomtype at,int *n,int typelist[],
+ int thistype,
+ t_param param[],int ftype)
+{
+ int i,nn,nrfp,j,k,ntype,tli;
+ bool bFound=FALSE;
+
+ nn = *n;
+ nrfp = NRFP(ftype);
+ ntype = get_atomtype_ntypes(at);
+
+ for(i=0; (i<nn); i++)
+ {
+ if (typelist[i] == thistype)
+ {
+ /* This type number has already been added */
+ break;
+ }
+
+ /* Otherwise, check if the parameters are identical to any previously added type */
+
+ bFound=TRUE;
+ for(j=0;j<ntype && bFound;j++)
+ {
+ /* Check nonbonded parameters */
+ for(k=0;k<nrfp && bFound;k++)
+ {
+ bFound=(param[ntype*typelist[i]+j].c[k]==param[ntype*thistype+j].c[k]);
+ }
+
+ /* Check radius, volume, surftens */
+ tli = typelist[i];
+ bFound = bFound &&
+ (get_atomtype_radius(tli,at) == get_atomtype_radius(thistype,at)) &&
+ (get_atomtype_vol(tli,at) == get_atomtype_vol(thistype,at)) &&
+ (get_atomtype_surftens(tli,at) == get_atomtype_surftens(thistype,at)) &&
+ (get_atomtype_atomnumber(tli,at) == get_atomtype_atomnumber(thistype,at));
+ }
+ if (bFound)
+ {
+ break;
+ }
+ }
+
+ if (i == nn) {
+ if (debug)
+ fprintf(debug,"Renumbering atomtype %d to %d\n",thistype,nn);
+ if (nn == ntype)
+ gmx_fatal(FARGS,"Atomtype horror n = %d, %s, %d",nn,__FILE__,__LINE__);
+ typelist[nn]=thistype;
+ nn++;
+ }
+ *n = nn;
+
+ return i;
+}
+
+void renum_atype(t_params plist[],t_topology *top,
+ int *wall_atomtype,
+ t_atomtype at,bool bVerbose)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) at;
+
+ int i,j,k,l,mi,mj,nat,nrfp,ftype,ntype;
+ t_param *nbsnew;
+ int *typelist;
+ real *new_radius;
+ real *new_vol;
+ real *new_surftens;
+ int *new_atomnumber;
+
+ ntype = get_atomtype_ntypes(at);
+ snew(typelist,ntype);
+
+ if (bVerbose)
+ fprintf(stderr,"renumbering atomtypes...\n");
+
+ /* Since the bonded interactions have been assigned now,
+ * we want to reduce the number of atom types by merging
+ * ones with identical nonbonded interactions, in addition
+ * to removing unused ones.
+ *
+ * With Generalized-Born electrostatics, or implicit solvent
+ * we also check that the atomtype radius, effective_volume
+ * and surface tension match.
+ *
+ * With QM/MM we also check that the atom numbers match
+ */
+
+ /* Get nonbonded interaction type */
+ if (plist[F_LJ].nr > 0)
+ ftype=F_LJ;
+ else
+ ftype=F_BHAM;
+
+ /* Renumber atomtypes by first making a list of which ones are actually used.
+ * We provide the list of nonbonded parameters so search_atomtypes
+ * can determine if two types should be merged.
+ */
+ nat=0;
+ for(i=0; (i<top->atoms.nr); i++) {
+ top->atoms.atom[i].type=
+ search_atomtypes(at,&nat,typelist,top->atoms.atom[i].type,
+ plist[ftype].param,ftype);
+ top->atoms.atom[i].typeB=
+ search_atomtypes(at,&nat,typelist,top->atoms.atom[i].typeB,
+ plist[ftype].param,ftype);
+ }
+
+ for(i=0; i<2; i++) {
+ if (wall_atomtype[i] >= 0)
+ wall_atomtype[i] = search_atomtypes(at,&nat,typelist,wall_atomtype[i],
+ plist[ftype].param,ftype);
+ }
+
+ snew(new_radius,nat);
+ snew(new_vol,nat);
+ snew(new_surftens,nat);
+ snew(new_atomnumber,nat);
+
+ /* We now have a list of unique atomtypes in typelist */
+
+ if (debug)
+ pr_ivec(debug,0,"typelist",typelist,nat,TRUE);
+
+ /* Renumber nlist */
+ nbsnew = NULL;
+ snew(nbsnew,plist[ftype].nr);
+
+ nrfp = NRFP(ftype);
+
+ for(i=k=0; (i<nat); i++)
+ {
+ mi=typelist[i];
+ for(j=0; (j<nat); j++,k++)
+ {
+ mj=typelist[j];
+ for(l=0; (l<nrfp); l++)
+ {
+ nbsnew[k].c[l]=plist[ftype].param[ntype*mi+mj].c[l];
+ }
+ }
+ new_radius[i] = get_atomtype_radius(mi,at);
+ new_vol[i] = get_atomtype_vol(mi,at);
+ new_surftens[i] = get_atomtype_surftens(mi,at);
+ new_atomnumber[i] = get_atomtype_atomnumber(mi,at);
+ }
+
+ for(i=0; (i<nat*nat); i++) {
+ for(l=0; (l<nrfp); l++)
+ plist[ftype].param[i].c[l]=nbsnew[i].c[l];
+ }
+ plist[ftype].nr=i;
+ top->idef.atnr = nat;
+
+ sfree(ga->radius);
+ sfree(ga->vol);
+ sfree(ga->surftens);
+ sfree(ga->atomnumber);
+
+ ga->radius = new_radius;
+ ga->vol = new_vol;
+ ga->surftens = new_surftens;
+ ga->atomnumber = new_atomnumber;
+
+ ga->nr=nat;
+
+ sfree(nbsnew);
+ sfree(typelist);
+}
+
+void copy_atomtype_atomtypes(t_atomtype atype,t_atomtypes *atomtypes)
+{
+ gpp_atomtype *ga = (gpp_atomtype *) atype;
+ int i,ntype;
+
+ /* Copy the atomtype data to the topology atomtype list */
+ ntype = get_atomtype_ntypes(atype);
+ atomtypes->nr=ntype;
+ snew(atomtypes->radius,ntype);
+ snew(atomtypes->vol,ntype);
+ snew(atomtypes->surftens,ntype);
+ snew(atomtypes->atomnumber,ntype);
+
+
+ for(i=0; i<ntype; i++) {
+ atomtypes->radius[i] = ga->radius[i];
+ atomtypes->vol[i] = ga->vol[i];
+ atomtypes->surftens[i] = ga->surftens[i];
+ atomtypes->atomnumber[i] = ga->atomnumber[i];
+ }
+}
+
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gromacs Runs On Most of All Computer Systems
+ */
+
+#ifndef _gpp_atomtype_h
+#define _gpp_atomtype_h
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include "typedefs.h"
+#include "macros.h"
+
+typedef struct gpp_atomtype *t_atomtype;
+
+extern int get_atomtype_type(char *str,t_atomtype at);
+/* Return atomtype corresponding to case-insensitive str
+ or NOTSET if not found */
+
+extern int get_atomtype_ntypes(t_atomtype at);
+/* Return number of atomtypes */
+
+extern char *get_atomtype_name(int nt, t_atomtype at);
+/* Return name corresponding to atomtype nt, or NULL if not found */
+
+extern real get_atomtype_massA(int nt,t_atomtype at);
+extern real get_atomtype_massB(int nt,t_atomtype at);
+extern real get_atomtype_qA(int nt,t_atomtype at);
+extern real get_atomtype_qB(int nt,t_atomtype at);
+extern real get_atomtype_radius(int nt,t_atomtype at);
+extern real get_atomtype_vol(int nt,t_atomtype at);
+extern real get_atomtype_surftens(int nt,t_atomtype at);
+extern int get_atomtype_ptype(int nt,t_atomtype at);
+extern int get_atomtype_batype(int nt,t_atomtype at);
+extern int get_atomtype_atomnumber(int nt,t_atomtype at);
+
+/* Return the above variable for atomtype nt, or NOTSET if not found */
+
+extern real get_atomtype_nbparam(int nt,int param,t_atomtype at);
+/* Similar to the previous but returns the paramth parameter or NOTSET */
+
+extern t_atomtype init_atomtype(void);
+/* Return a new atomtype structure */
+
+extern void done_atomtype(t_atomtype *at);
+/* Free the memory in the structure */
+
+extern int set_atomtype(int nt,t_atomtype at,t_symtab *tab,
+ t_atom *a,char *name,t_param *nb,
+ int bondatomtype,
+ real radius,real vol,real surftens,int atomnumber);
+/* Set the values of an existing atom type nt. Returns nt on success or
+ NOTSET on error. */
+
+extern int add_atomtype(t_atomtype at,t_symtab *tab,
+ t_atom *a,char *name,t_param *nb,
+ int bondatomtype,
+ real radius,real vol,real surftens,real atomnumber);
+/* Add a complete new atom type to an existing atomtype structure. Returns
+ the number of the atom type. */
+
+extern void print_at (FILE * out, t_atomtype at);
+/* Print an atomtype record to a text file */
+
+extern void renum_atype(t_params plist[],t_topology *top,
+ int *wall_atomtype,
+ t_atomtype at,bool bVerbose);
+
+extern void copy_atomtype_atomtypes(t_atomtype atype,t_atomtypes *atypes);
+/* Copy from one structure to another */
+#endif /* _gpp_atomtype_h */
+
+
+
+
+
+
+
+
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "smalloc.h"
+#include "sysstuff.h"
+#include "macros.h"
+#include "symtab.h"
+#include "string2.h"
+#include "gpp_bond_atomtype.h"
+
+typedef struct {
+ int nr; /* The number of atomtypes */
+ char ***atomname; /* Names of the atomtypes */
+} gpp_bond_atomtype;
+
+int get_bond_atomtype_type(char *str,t_bond_atomtype at)
+{
+ gpp_bond_atomtype *ga = (gpp_bond_atomtype *) at;
+
+ int i;
+
+ for (i=0; (i<ga->nr); i++)
+ if (strcasecmp(str,*(ga->atomname[i])) == 0)
+ return i;
+
+ return NOTSET;
+}
+
+char *get_bond_atomtype_name(int nt, t_bond_atomtype at)
+{
+ gpp_bond_atomtype *ga = (gpp_bond_atomtype *) at;
+
+ if ((nt < 0) || (nt >= ga->nr))
+ return NULL;
+
+ return *(ga->atomname[nt]);
+}
+
+t_bond_atomtype init_bond_atomtype(void)
+{
+ gpp_bond_atomtype *ga;
+
+ snew(ga,1);
+
+ return (t_bond_atomtype ) ga;
+}
+
+void add_bond_atomtype(t_bond_atomtype at,t_symtab *tab,
+ char *name)
+{
+ gpp_bond_atomtype *ga = (gpp_bond_atomtype *) at;
+
+ ga->nr++;
+ srenew(ga->atomname,ga->nr);
+ ga->atomname[ga->nr-1] = put_symtab(tab,name);
+}
+
+void done_bond_atomtype(t_bond_atomtype *at)
+{
+ gpp_bond_atomtype *ga = (gpp_bond_atomtype *) *at;
+
+ sfree(ga->atomname);
+ ga->nr = 0;
+ sfree(ga);
+
+ *at = NULL;
+}
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gromacs Runs On Most of All Computer Systems
+ */
+
+#ifndef _gpp_bondatomtype_h
+#define _gpp_bondatomtype_h
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include "typedefs.h"
+#include "macros.h"
+
+typedef struct gpp_bondatomtype *t_bond_atomtype;
+
+extern int get_bond_atomtype_type(char *str,t_bond_atomtype at);
+/* Return atomtype corresponding to case-insensitive str
+ or NOTSET if not found */
+
+extern char *get_bond_atomtype_name(int nt, t_bond_atomtype at);
+/* Return name corresponding to atomtype nt, or NULL if not found */
+
+extern t_bond_atomtype init_bond_atomtype(void);
+/* Return a new atomtype structure */
+
+extern void done_bond_atomtype(t_bond_atomtype *at);
+/* Free the memory in the structure */
+
+extern void add_bond_atomtype(t_bond_atomtype at,t_symtab *tab,
+ char *name);
+/* Add a complete new atom type to an existing atomtype structure */
+
+#endif /* _gpp_bond_atomtype_h */
+
+
+
+
+
+
+
+
--- /dev/null
+/*
+ * $Id$
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gromacs Runs On Most of All Computer Systems
+ */
+
+#ifndef _tomorse_h
+#define _tomorse_h
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include "typedefs.h"
+#include "macros.h"
+
+extern void convert_harmonics(int nrmols,t_molinfo mols[],t_atomtype atype);
+
+#endif /* _grompp_h */
+
+
+
+
+
+
+
+
#include "enxio.h"
#include "perf_est.h"
#include "compute_io.h"
+#include "gpp_atomtype.h"
+#include "gpp_tomorse.h"
static int rm_interactions(int ifunc,int nrmols,t_molinfo mols[])
{
}
}
-static void check_pairs(int nrmols,t_molinfo mi[],int ntype,t_param *nb)
-{
- int i,j,jj,k,ai,aj,taiA,tajA,taiB,tajB,bLJ;
- t_params *p;
-
- for(i=0; (i<nrmols); i++) {
- p = &(mi[i].plist[F_LJ14]);
- for(j=jj=0; (j<p->nr); j++) {
- /* Extract atom types and hence the index in the nbf matrix */
- ai = p->param[j].a[0];
- aj = p->param[j].a[1];
- taiA = mi[i].atoms.atom[ai].type;
- tajA = mi[i].atoms.atom[aj].type;
- taiB = mi[i].atoms.atom[ai].typeB;
- tajB = mi[i].atoms.atom[aj].typeB;
-
- bLJ = FALSE;
- for(k=0; (k<MAXFORCEPARAM); k++)
- bLJ = bLJ || (((nb[taiA].c[k]*nb[tajA].c[k]) != 0) ||
- ((nb[taiB].c[k]*nb[tajB].c[k]) != 0));
- if (bLJ) {
- cp_param(&(p->param[jj]),&(p->param[j]));
- jj++;
- }
- else if (debug) {
- fprintf(debug,"Removed 1-4 interaction between atoms %d and %d (within mol %s)\n",
- ai+1,aj+1,*(mi[i].name));
- }
- }
- fprintf(stderr,"Removed %d 1-4 interactions for molecule %s\n",
- p->nr-jj,*(mi[i].name));
- p->nr = jj;
- }
-}
-
static void check_vel(t_atoms *atoms,rvec v[])
{
int i;
new_status(char *topfile,char *topppfile,char *confin,
t_gromppopts *opts,t_inputrec *ir,bool bZero,
bool bGenVel,bool bVerbose,t_state *state,
- t_atomtype *atype,t_topology *sys,
+ t_atomtype atype,t_topology *sys,
t_molinfo *msys,t_params plist[],int *comb,real *reppow,
bool bEnsemble,bool bMorse,
- bool bCheckPairs,int *nerror)
+ int *nerror)
{
t_molinfo *molinfo=NULL;
t_simsystem *Sims=NULL;
msys->name=do_top(bVerbose,topfile,topppfile,opts,bZero,&(sys->symtab),
plist,comb,reppow,atype,&nrmols,&molinfo,ir,&Nsim,&Sims);
- if (bCheckPairs)
- check_pairs(nrmols,molinfo,atype->nr,atype->nb);
-
if (bMorse)
convert_harmonics(nrmols,molinfo,atype);
}
}
-static int search_atomtypes(t_atomtype *at,int *n,int typelist[],int thistype,
- t_param param[],int ftype)
-{
- int i,nn,nrfp,j,k,found;
-
- nn = *n;
- nrfp = NRFP(ftype);
-
- for(i=0; (i<nn); i++)
- {
- if (typelist[i] == thistype)
- {
- /* This type number has already been added */
- break;
- }
-
- /* Otherwise, check if the parameters are identical to any previously added type */
-
- found=1;
- for(j=0;j<at->nr && found;j++)
- {
- /* Check nonbonded parameters */
- for(k=0;k<nrfp && found;k++)
- {
- found=(param[at->nr*typelist[i]+j].c[k]==param[at->nr*thistype+j].c[k]);
- }
-
- /* Check radius, volume, surftens */
- found = found &&
- ((at->radius[typelist[i]] == at->radius[thistype]) &&
- (at->vol[typelist[i]] == at->vol[thistype]) &&
- (at->surftens[typelist[i]] == at->surftens[thistype]) &&
- (at->atomnumber[typelist[i]] == at->atomnumber[thistype]));
- }
- if (found)
- {
- break;
- }
- }
-
- if (i == nn) {
- if (debug)
- fprintf(debug,"Renumbering atomtype %d to %d\n",thistype,nn);
- if (nn == at->nr)
- gmx_fatal(FARGS,"Atomtype horror n = %d, %s, %d",nn,__FILE__,__LINE__);
- typelist[nn]=thistype;
- nn++;
- }
- *n = nn;
-
- return i;
-}
-
-static int renum_atype(t_params plist[],t_topology *top,int *wall_atomtype,
- t_atomtype *at,bool bVerbose)
-{
- int i,j,k,l,mi,mj,nat,nrfp,ftype;
- t_param *nbsnew;
- int *typelist;
- real *new_radius;
- real *new_vol;
- real *new_surftens;
- int *new_atomnumber;
-
- snew(typelist,at->nr);
-
- if (bVerbose)
- fprintf(stderr,"renumbering atomtypes...\n");
-
- /* Since the bonded interactions have been assigned now,
- * we want to reduce the number of atom types by merging
- * ones with identical nonbonded interactions, in addition
- * to removing unused ones.
- *
- * With Generalized-Born electrostatics, or implicit solvent
- * we also check that the atomtype radius, effective_volume
- * and surface tension match.
- *
- * With QM/MM we also check that the atom numbers match
- */
-
- /* Get nonbonded interaction type */
- if (plist[F_LJ].nr > 0)
- ftype=F_LJ;
- else
- ftype=F_BHAM;
-
- /* Renumber atomtypes by first making a list of which ones are actually used.
- * We provide the list of nonbonded parameters so search_atomtypes
- * can determine if two types should be merged.
- */
- nat=0;
- for(i=0; (i<top->atoms.nr); i++) {
- top->atoms.atom[i].type=
- search_atomtypes(at,&nat,typelist,top->atoms.atom[i].type,
- plist[ftype].param,ftype);
- top->atoms.atom[i].typeB=
- search_atomtypes(at,&nat,typelist,top->atoms.atom[i].typeB,
- plist[ftype].param,ftype);
- }
-
- for(i=0; i<2; i++) {
- if (wall_atomtype[i] >= 0)
- wall_atomtype[i] = search_atomtypes(at,&nat,typelist,wall_atomtype[i],
- plist[ftype].param,ftype);
- }
-
- snew(new_radius,nat);
- snew(new_vol,nat);
- snew(new_surftens,nat);
- snew(new_atomnumber,nat);
-
- /* We now have a list of unique atomtypes in typelist */
-
- if (debug)
- pr_ivec(debug,0,"typelist",typelist,nat,TRUE);
-
- /* Renumber nlist */
- nbsnew = NULL;
- snew(nbsnew,plist[ftype].nr);
-
- nrfp = NRFP(ftype);
-
- for(i=k=0; (i<nat); i++)
- {
- mi=typelist[i];
- for(j=0; (j<nat); j++,k++)
- {
- mj=typelist[j];
- for(l=0; (l<nrfp); l++)
- {
- nbsnew[k].c[l]=plist[ftype].param[at->nr*mi+mj].c[l];
- }
- }
- new_radius[i] = at->radius[mi];
- new_vol[i] = at->vol[mi];
- new_surftens[i] = at->surftens[mi];
- new_atomnumber[i] = at->atomnumber[mi];
- }
-
- for(i=0; (i<nat*nat); i++) {
- for(l=0; (l<nrfp); l++)
- plist[ftype].param[i].c[l]=nbsnew[i].c[l];
- }
- plist[ftype].nr=i;
-
- sfree(at->radius);
- sfree(at->vol);
- sfree(at->surftens);
- sfree(at->atomnumber);
-
- at->radius = new_radius;
- at->vol = new_vol;
- at->surftens = new_surftens;
- at->atomnumber = new_atomnumber;
-
- at->nr=nat;
-
- sfree(nbsnew);
- sfree(typelist);
-
- return nat;
-}
-
-static void set_wall_atomtype(t_atomtype *at,t_gromppopts *opts,
+static void set_wall_atomtype(t_atomtype at,t_gromppopts *opts,
t_inputrec *ir)
{
int i;
if (ir->nwall > 0)
fprintf(stderr,"Searching the wall atom type(s)\n");
for(i=0; i<ir->nwall; i++)
- ir->wall_atomtype[i] = at2type(opts->wall_atomtype[i],at);
+ ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i],at);
}
static int count_constraints(t_params plist[])
t_gromppopts *opts;
t_topology *sys;
t_molinfo msys;
- t_atomtype *atype;
+ t_atomtype atype;
t_inputrec *ir;
int natoms,nvsite,comb;
t_params *plist;
matrix box;
real max_spacing,reppow;
char fn[STRLEN],fnB[STRLEN],*mdparin;
- int nerror;
+ int nerror,ntype;
bool bNeedVel,bGenVel;
bool have_radius,have_vol,have_surftens;
bool have_atomnumber;
/* Command line options */
static bool bVerbose=TRUE,bRenum=TRUE;
- static bool bRmVSBds=TRUE,bCheckPairs=FALSE,bZero=FALSE;
+ static bool bRmVSBds=TRUE,bZero=FALSE;
static int i,maxwarn=0;
static real fr_time=-1;
t_pargs pa[] = {
"Remove constant bonded interactions with virtual sites" },
{ "-maxwarn", FALSE, etINT, {&maxwarn},
"Number of allowed warnings during input processing" },
- { "-check14", FALSE, etBOOL, {&bCheckPairs},
- "Remove 1-4 interactions without Van der Waals" },
{ "-zero", FALSE, etBOOL, {&bZero},
"Set parameters for bonded interactions without defaults to zero instead of generating an error" },
{ "-renum", FALSE, etBOOL, {&bRenum},
snew(plist,F_NRE);
init_plist(plist);
snew(sys,1);
- snew(atype,1);
+ atype = init_atomtype();
if (debug)
pr_symtab(debug,0,"Just opened",&sys->symtab);
opts,ir,bZero,bGenVel,bVerbose,&state,
atype,sys,&msys,plist,&comb,&reppow,
(opts->eDisre==edrEnsemble),opts->bMorse,
- bCheckPairs,&nerror);
+ &nerror);
if (debug)
pr_symtab(debug,0,"After new_status",&sys->symtab);
/* If we are doing GBSA, check that we got the parameters we need */
have_radius=have_vol=have_surftens=TRUE;
- for(i=0;i<atype->nr;i++) {
- have_radius=have_radius && (atype->radius[i]>0);
- have_vol=have_vol && (atype->vol[i]>0);
- have_surftens=have_surftens && (atype->surftens[i]>=0);
+ ntype = get_atomtype_ntypes(atype);
+ for(i=0; (i<ntype); i++) {
+ have_radius = have_radius && (get_atomtype_radius(i,atype) > 0);
+ have_vol = have_vol && (get_atomtype_vol(i,atype) > 0);
+ have_surftens = have_surftens && (get_atomtype_surftens(i,atype) >= 0);
}
- if(!have_radius && ir->coulombtype==eelGB) {
+ if (!have_radius && ir->coulombtype==eelGB) {
fprintf(stderr,"Can't do GB electrostatics; the forcefield is missing values for\n"
"atomtype radii, or they might be zero.");
nerror++;
/* If we are doing QM/MM, check that we got the atom numbers */
have_atomnumber = TRUE;
- for (i=0;i<atype->nr;i++) {
- have_atomnumber = have_atomnumber && (atype->atomnumber[i]>=0);
+ for (i=0; i<ntype; i++) {
+ have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i,atype) >= 0);
}
if (!have_atomnumber && ir->bQMMM)
{
clean_vsite_bondeds(msys.plist,sys->atoms.nr,bRmVSBds);
set_wall_atomtype(atype,opts,ir);
- if (bRenum)
- atype->nr = renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
+ if (bRenum) {
+ renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
+ ntype = get_atomtype_ntypes(atype);
+ }
/* Copy the atomtype data to the topology atomtype list */
- sys->atomtypes.nr=atype->nr;
- snew(sys->atomtypes.radius,atype->nr);
- snew(sys->atomtypes.vol,atype->nr);
- snew(sys->atomtypes.surftens,atype->nr);
- snew(sys->atomtypes.atomnumber,atype->nr);
-
-
- for(i=0;i<atype->nr;i++) {
- sys->atomtypes.radius[i]=atype->radius[i];
- sys->atomtypes.vol[i]=atype->vol[i];
- sys->atomtypes.surftens[i]=atype->surftens[i];
- sys->atomtypes.atomnumber[i] = atype->atomnumber[i];
- }
+ copy_atomtype_atomtypes(atype,&(sys->atomtypes));
if (debug)
pr_symtab(debug,0,"After renum_atype",&sys->symtab);
if (bVerbose)
fprintf(stderr,"converting bonded parameters...\n");
- convert_params(atype->nr, plist, msys.plist, comb, reppow, &sys->idef);
+ convert_params(ntype, plist, msys.plist, comb, reppow, &sys->idef);
if (debug)
pr_symtab(debug,0,"After convert_params",&sys->symtab);
#include "typedefs.h"
#include "pdbio.h"
#include "grompp.h"
+#include "gpp_atomtype.h"
/* Used for reading .rtp/.tdb */
/* ebtsBONDS must be the first, new types can be added to the end */
int nah;
t_symtab tab;
int *rN, *rC;
- t_atomtype *atype;
+ t_atomtype atype;
/* protonated topology: */
t_atoms *patoms;
/* unprotonated topology: */
}
int nm2type(x2top_nm2t nm2t,t_symtab *tab,t_atoms *atoms,
- t_atomtype *atype,int *nbonds,t_params *bonds)
+ t_atomtype atype,int *nbonds,t_params *bonds)
{
x2top_nm2type *nm2type = (x2top_nm2type *) nm2t;
- int cur = 0;
-#define prev (1-cur)
- int i,j,k,m,n,nresolved,nb,maxbond,ai,aj,best,im,nqual[2][ematchNR];
- int *bbb,*n_mask,*m_mask,**match,**quality;
- char *aname_i,*aname_m,*aname_n,*type;
- double qq,mm;
-
+ int cur = 0;
+#define prev (1-cur)
+ int i,j,k,m,n,nresolved,nb,maxbond;
+ int ai,aj,best,im,nqual[2][ematchNR];
+ int *bbb,*n_mask,*m_mask,**match,**quality;
+ char *aname_i,*aname_m,*aname_n,*type;
+ double alpha,mm;
+ t_atom *atom;
+ t_param *param;
+
+ snew(atom,1);
+ snew(param,1);
maxbond = 0;
for(i=0; (i<atoms->nr); i++)
maxbond = max(maxbond,nbonds[i]);
}
}
if (best != -1) {
- qq = nm2type->nm2t[best].alpha;
- mm = nm2type->nm2t[best].m;
- type = nm2type->nm2t[best].type;
-
- for(k=0; (k<atype->nr); k++) {
- if (strcasecmp(*atype->atomname[k],type) == 0)
- break;
+ alpha = nm2type->nm2t[best].alpha;
+ mm = nm2type->nm2t[best].m;
+ type = nm2type->nm2t[best].type;
+
+ if ((k = get_atomtype_type(type,atype)) == NOTSET) {
+ atom->type = atom->typeB = get_atomtype_ntypes(atype);
+ atom->qB = alpha;
+ atom->m = atom->mB = mm;
+ add_atomtype(atype,tab,atom,type,param,
+ atom->type,0,0,0,0);
+ }
+ else {
+ atoms->atom[i].type = k;
+ atoms->atom[i].typeB = k;
+ atoms->atom[i].q = 0;
+ atoms->atom[i].qB = alpha;
+ atoms->atom[i].m = mm;
+ atoms->atom[i].mB = mm;
}
- if (k == atype->nr) {
- /* New atomtype */
- atype->nr++;
- srenew(atype->atomname,atype->nr);
- srenew(atype->atom,atype->nr);
- srenew(atype->bondatomtype,atype->nr);
- atype->atomname[k] = put_symtab(tab,type);
- atype->bondatomtype[k] = k; /* Set bond_atomtype identical to atomtype */
- atype->atom[k].type = k;
- atype->atom[k].typeB = k;
- atype->atom[k].q = 0;
- atype->atom[k].qB = qq;
- atype->atom[k].m = mm;
- atype->atom[k].mB = mm;
- }
- atoms->atom[i].type = k;
- atoms->atom[i].typeB = k;
- atoms->atom[i].q = 0;
- atoms->atom[i].qB = qq;
- atoms->atom[i].m = mm;
- atoms->atom[i].mB = mm;
nresolved++;
}
else {
}
sfree(bbb);
sfree(n_mask);
+ for(i=0; (i<maxbond); i++)
+ sfree(match[i]);
+ sfree(match);
sfree(m_mask);
-
+ sfree(atom);
+ sfree(param);
+
return nresolved;
}
t_restp *restp;
t_hackblock *ah;
t_symtab symtab;
- t_atomtype *atype;
+ t_atomtype atype;
t_aa_names *aan;
char fn[256],*top_fn,itp_fn[STRLEN],posre_fn[STRLEN],buf_fn[STRLEN];
char molname[STRLEN],title[STRLEN],resname[STRLEN],quote[STRLEN];
}
-static int name2type(t_atoms *at, int **cgnr, t_atomtype *atype,
+static int name2type(t_atoms *at, int **cgnr, t_atomtype atype,
t_restp restp[])
{
int i,j,prevresnr,resnr,i0,prevcg,cg,curcg;
j=search_jtype(&restp[resnr],name,bNterm);
at->atom[i].type = restp[resnr].atom[j].type;
at->atom[i].q = restp[resnr].atom[j].q;
- at->atom[i].m = atype->atom[restp[resnr].atom[j].type].m;
+ at->atom[i].m = get_atomtype_massA(restp[resnr].atom[j].type,
+ atype);
cg = restp[resnr].cgnr[j];
if ( (cg != prevcg) || (resnr != prevresnr) )
curcg++;
void write_top(FILE *out, char *pr,char *molname,
t_atoms *at,int bts[],t_params plist[],t_excls excls[],
- t_atomtype *atype,int *cgnr, int nrexcl)
+ t_atomtype atype,int *cgnr, int nrexcl)
/* NOTE: nrexcl is not the size of *excl! */
{
if (at && atype && cgnr) {
}
void pdb2top(FILE *top_file, char *posre_fn, char *molname,
- t_atoms *atoms, rvec **x, t_atomtype *atype, t_symtab *tab,
+ t_atoms *atoms, rvec **x, t_atomtype atype, t_symtab *tab,
int bts[], int nrtp, t_restp rtp[],
int nterpairs,t_hackblock **ntdb, t_hackblock **ctdb,
int *rn, int *rc, bool bMissing,
extern void pdb2top(FILE *top_file, char *posre_fn, char *molname,
t_atoms *atoms,rvec **x,
- t_atomtype *atype,t_symtab *tab,
+ t_atomtype atype,t_symtab *tab,
int bts[],
int nrtp, t_restp rtp[],
int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
extern void write_top(FILE *out, char *pr,char *molname,
t_atoms *at,int bts[],t_params plist[],t_excls excls[],
- t_atomtype *atype,int *cgnr, int nrexcl);
+ t_atomtype atype,int *cgnr, int nrexcl);
/* NOTE: nrexcl is not the size of *excl! */
#endif /* _pdb2top_h */
#include "resall.h"
#include "pgutil.h"
-t_atomtype *read_atype(char *adb,t_symtab *tab)
+t_atomtype read_atype(char *adb,t_symtab *tab)
{
-#define MAXAT 5000
FILE *in;
char aadb[STRLEN];
char buf[STRLEN],name[STRLEN];
double m;
- int nratt;
- t_atomtype *at;
-
- sprintf(aadb,"%s.atp",adb);
- in=libopen(aadb);
+ int nratt=0;
+ t_atomtype at;
+ t_atom *a;
+ t_param *nb;
- snew(at,1);
- snew(at->atom,MAXAT);
- snew(at->atomname,MAXAT);
+ sprintf(aadb,"%s.atp",adb);
+ in = libopen(aadb);
+ at = init_atomtype();
+ snew(a,1);
+ snew(nb,1);
- for(nratt=0; ; nratt++) {
- if (nratt >= MAXAT)
- gmx_fatal(FARGS,"nratt >= MAXAT(%d). Increase the latter",MAXAT);
- if (feof(in))
- break;
-
+ while(!feof(in)) {
/* Skip blank or comment-only lines */
do {
fgets2(buf,STRLEN,in);
- if(buf) {
+ if (buf) {
strip_comment(buf);
trim(buf);
}
} while (buf && strlen(buf)==0);
-
- if(buf==NULL)
- break;
- if (sscanf(buf,"%s%lf",name,&m) != 2)
- break;
- set_at(&(at->atom[nratt]),m,0.0,nratt,0);
- at->atomname[nratt]=put_symtab(tab,name);
- fprintf(stderr,"\rAtomtype %d",nratt+1);
+ if ((buf != NULL) && (sscanf(buf,"%s%lf",name,&m) == 2)) {
+ a->m = m;
+ add_atomtype(at,tab,a,name,nb,0,0,0,0,0);
+ fprintf(stderr,"\rAtomtype %d",nratt+1);
+ }
}
fclose(in);
fprintf(stderr,"\n");
- at->nr=nratt;
return at;
}
-static void print_resatoms(FILE *out,t_atomtype *atype,t_restp *rtp)
+static void print_resatoms(FILE *out,t_atomtype atype,t_restp *rtp)
{
int j,tp;
+ char *tpnm;
/* fprintf(out,"%5s\n",rtp->resname);
fprintf(out,"%5d\n",rtp->natom); */
fprintf(out," [ atoms ]\n");
for(j=0; (j<rtp->natom); j++) {
- tp=rtp->atom[j].type;
- if ((tp < 0) || (tp >= atype->nr))
- gmx_fatal(FARGS,"tp (%d) out of range (0 .. %d)",tp,atype->nr);
+ tp = rtp->atom[j].type;
+ tpnm = get_atomtype_name(tp,atype);
+ if (tpnm == NULL)
+ gmx_fatal(FARGS,"Incorrect atomtype (%d)",tp);
fprintf(out,"%6s%6s%8.3f%6d\n",
- *(rtp->atomname[j]),*(atype->atomname[tp]),
- rtp->atom[j].q,rtp->cgnr[j]);
+ *(rtp->atomname[j]),tpnm,rtp->atom[j].q,rtp->cgnr[j]);
}
}
static bool read_atoms(FILE *in,char *line,
- t_restp *r0,t_symtab *tab,t_atomtype *atype)
+ t_restp *r0,t_symtab *tab,t_atomtype atype)
{
int i,j,cg,maxentries;
char buf[256],buf1[256];
r0->atomname[i]=put_symtab(tab,buf);
r0->atom[i].q=q;
r0->cgnr[i]=cg;
- for(j=0; (j<atype->nr); j++)
- if (strcasecmp(buf1,*(atype->atomname[j])) == 0)
- break;
- if (j == atype->nr)
+ j = get_atomtype_type(buf1,atype);
+ if (j == NOTSET)
gmx_fatal(FARGS,"Atom type %s (residue %s) not found in atomtype "
"database",buf1,r0->resname);
r0->atom[i].type=j;
- r0->atom[i].m=atype->atom[j].m;
+ r0->atom[i].m=get_atomtype_massA(j,atype);
i++;
}
r0->natom=i;
}
int read_resall(char *ff, int bts[], t_restp **rtp,
- t_atomtype *atype, t_symtab *tab, bool *bAlldih, int *nrexcl,
+ t_atomtype atype, t_symtab *tab, bool *bAlldih, int *nrexcl,
bool *HH14, bool *bRemoveDih)
{
FILE *in;
}
void print_resall(FILE *out, int bts[], int nrtp, t_restp rtp[],
- t_atomtype *atype, bool bAlldih, int nrexcl,
+ t_atomtype atype, bool bAlldih, int nrexcl,
bool HH14, bool bRemoveDih)
{
int i,bt;
#include "typedefs.h"
#include "hackblock.h"
+#include "gpp_atomtype.h"
#include "grompp.h"
extern t_restp *search_rtp(char *key,int nrtp,t_restp rtp[]);
/* Search for an entry in the rtp database */
-extern t_atomtype *read_atype(char *adb,t_symtab *tab);
+extern t_atomtype read_atype(char *adb,t_symtab *tab);
/* read atom type database */
extern int read_resall(char *resdb, int bts[], t_restp **rtp,
- t_atomtype *atype, t_symtab *tab, bool *bAlldih, int *nrexcl,
+ t_atomtype atype, t_symtab *tab, bool *bAlldih, int *nrexcl,
bool *HH14, bool *bRemoveDih);
/* read rtp database */
extern void print_resall(FILE *out, int bts[], int nrtp, t_restp rtp[],
- t_atomtype *atype, bool bAlldih, int nrexcl,
+ t_atomtype atype, bool bAlldih, int nrexcl,
bool HH14, bool remove_dih);
/* write rtp database */
#define FATAL() gmx_fatal(FARGS,"Reading Termini Database: not enough items on line\n%s",line)
-static void read_atom(char *line, t_atom *a, t_atomtype *atype, int *cgnr)
+static void read_atom(char *line, t_atom *a, t_atomtype atype, int *cgnr)
{
int i, n;
char type[12];
gmx_fatal(FARGS,"Reading Termini Database: expected %d items of atom data in stead of %d on line\n%s", 3, i, line);
a->m=m;
a->q=q;
- a->type=at2type(type,atype);
+ a->type=get_atomtype_type(type,atype);
if ( sscanf(line+n,"%d", cgnr) != 1 )
*cgnr = NOTSET;
}
-static void print_atom(FILE *out,t_atom *a,t_atomtype *atype,char *newnm)
+static void print_atom(FILE *out,t_atom *a,t_atomtype atype,char *newnm)
{
fprintf(out,"\t%s\t%g\t%g\n",
- type2nm(a->type,atype),a->m,a->q);
+ get_atomtype_name(a->type,atype),a->m,a->q);
}
-static void print_ter_db(char *ff,char C,int nb,t_hackblock tb[],t_atomtype *atype)
+static void print_ter_db(char *ff,char C,int nb,t_hackblock tb[],t_atomtype atype)
{
FILE *out;
int i,j,k,bt,nrepl,nadd,ndel;
fclose(out);
}
-int read_ter_db(char *FF,char ter,t_hackblock **tbptr,t_atomtype *atype)
+int read_ter_db(char *FF,char ter,t_hackblock **tbptr,t_atomtype atype)
{
FILE *in;
char inf[STRLEN],header[STRLEN],buf[STRLEN],line[STRLEN];
#include "grompp.h"
extern int read_ter_db(char *FF,char ter,
- t_hackblock **tbptr,t_atomtype *atype);
+ t_hackblock **tbptr,t_atomtype atype);
/* Read database for N&C terminal hacking */
extern t_hackblock **filter_ter(int nb,t_hackblock tb[],char *resname,int *nret);
#include "smalloc.h"
#include "toputil.h"
#include "gmx_fatal.h"
+#include "gpp_atomtype.h"
+#include "gpp_tomorse.h"
typedef struct {
char *ai,*aj;
}
}
-void convert_harmonics(int nrmols,t_molinfo mols[],t_atomtype *atype)
+void convert_harmonics(int nrmols,t_molinfo mols[],t_atomtype atype)
{
int n2m;
t_2morse *t2m;
for(j=0; (j<nrharm); j++) {
ni = mols[i].plist[bb].param[j].AI;
nj = mols[i].plist[bb].param[j].AJ;
- edis = search_e_diss(n2m,t2m,
- *atype->atomname[mols[i].atoms.atom[ni].type],
- *atype->atomname[mols[i].atoms.atom[nj].type]);
+ edis =
+ search_e_diss(n2m,t2m,
+ get_atomtype_name(mols[i].atoms.atom[ni].type,atype),
+ get_atomtype_name(mols[i].atoms.atom[nj].type,atype));
if (edis != 0) {
bRemoveHarm[j] = TRUE;
b0 = mols[i].plist[bb].param[j].c[0];
#include "topio.h"
#include "topshake.h"
#include "gmxcpp.h"
+#include "gpp_bond_atomtype.h"
#define CPPMARK '#' /* mark from cpp */
#define OPENDIR '[' /* starting sign for directive */
static char **read_topol(char *infile,char *outfile,
char *define,char *include,
t_symtab *symtab,
- t_atomtype *atype,
+ t_atomtype atype,
int *nrmols,
t_molinfo **molinfo,
t_params plist[],
real fudgeLJ=-1; /* Multiplication factor to generate 1-4 from LJ */
bool bReadDefaults,bReadMolType,bGenPairs,bWarn_copy_A_B;
double qt=0,qBt=0; /* total charge */
- t_bond_atomtype *batype;
+ t_bond_atomtype batype;
int lastcg=-1;
/* File handling variables */
int status,done;
bWarn_copy_A_B = bFEP;
- snew(batype,1);
- init_bond_atomtype(batype);
+ batype = init_bond_atomtype();
/* parse the actual file */
bReadDefaults = FALSE;
bGenPairs = FALSE;
&nbparam,bGenPairs ? &pair : NULL);
break;
-#define PUSHBT(nral) push_bt(d,plist,nral,batype->atomname,batype->nr,pline)
case d_bondtypes:
- PUSHBT(2);
+ push_bt(d,plist,2,NULL,batype,pline);
break;
case d_constrainttypes:
- PUSHBT(2);
+ push_bt(d,plist,2,NULL,batype,pline);
break;
case d_pairtypes:
if (bGenPairs)
push_nbt(d,pair,atype,pline,F_LJ14);
else
- push_bt(d,plist,2,atype->atomname,atype->nr,pline);
+ push_bt(d,plist,2,atype,NULL,pline);
break;
case d_angletypes:
- PUSHBT(3);
+ push_bt(d,plist,3,NULL,batype,pline);
break;
case d_dihedraltypes:
/* Special routine that can read both 2 and 4 atom dihedral definitions. */
- push_dihedraltype(d,plist,batype->atomname,batype->nr,pline);
+ push_dihedraltype(d,plist,batype,pline);
break;
-#undef PUSHBT
+
case d_nonbond_params:
push_nbt(d,nbparam,atype,pline,nb_funct);
break;
*/
case d_moleculetype: {
if (!bReadMolType) {
- ncombs = atype->nr*(atype->nr+1)/2;
+ int ntype = get_atomtype_ntypes(atype);
+ ncombs = (ntype*(ntype+1))/2;
generate_nbparams(comb,nb_funct,&(plist[nb_funct]),atype);
ncopy = copy_nbparams(nbparam,nb_funct,&(plist[nb_funct]),
- atype->nr);
+ ntype);
fprintf(stderr,"Generated %d of the %d non-bonded parameter combinations\n",ncombs-ncopy,ncombs);
- free_nbparam(nbparam,atype->nr);
+ free_nbparam(nbparam,ntype);
if (bGenPairs) {
gen_pairs(&(plist[nb_funct]),&(plist[F_LJ14]),fudgeLJ,comb,bVerbose);
ncopy = copy_nbparams(pair,nb_funct,&(plist[F_LJ14]),
- atype->nr);
+ ntype);
fprintf(stderr,"Generated %d of the %d 1-4 parameter combinations\n",ncombs-ncopy,ncombs);
- free_nbparam(pair,atype->nr);
+ free_nbparam(pair,ntype);
}
/* Copy GBSA parameters to atomtype array */
for(i=0; i<nmol; i++)
done_block2(&(block2[i]));
free(block2);
+
+ done_bond_atomtype(&batype);
+
*nrmols=nmol;
*nsim=Nsim;
t_params plist[],
int *combination_rule,
real *repulsion_power,
- t_atomtype *atype,
+ t_atomtype atype,
int *nrmols,
t_molinfo **molinfo,
t_inputrec *ir,
char *tmpfile;
char **title;
- init_atomtype(atype);
-
if (topppfile)
tmpfile = topppfile;
else
#include "typedefs.h"
#include "readir.h"
#include "grompp.h"
+#include "gpp_atomtype.h"
typedef struct {
int whichmol;
t_params plist[],
int *combination_rule,
real *repulsion_power,
- t_atomtype *atype,
+ t_atomtype atype,
int *nrmols,
t_molinfo **molinfo,
t_inputrec *ir,
#include "topdirs.h"
#include "symtab.h"
#include "gmx_fatal.h"
+#include "gpp_atomtype.h"
+#include "gpp_bond_atomtype.h"
-
-void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype *atype)
+void generate_nbparams(int comb,int ftype,t_params *plist,t_atomtype atype)
{
int i,j,k=-1,nf;
int nr,nrfp;
- real c,bi,bj;
+ real c,bi,bj,ci,cj,ci0,ci1,ci2,cj0,cj1,cj2;
char errbuf[256];
/* Lean mean shortcuts */
- nr = atype->nr;
+ nr = get_atomtype_ntypes(atype);
nrfp = NRFP(ftype);
snew(plist->param,nr*nr);
plist->nr=nr*nr;
switch (comb) {
case eCOMB_GEOMETRIC:
/* Gromos rules */
- for(i=k=0; (i<nr); i++)
- for(j=0; (j<nr); j++,k++)
+ for(i=k=0; (i<nr); i++) {
+ for(j=0; (j<nr); j++,k++) {
for(nf=0; (nf<nrfp); nf++) {
- c = sqrt(atype->nb[i].c[nf] * atype->nb[j].c[nf]);
+ ci = get_atomtype_nbparam(i,nf,atype);
+ cj = get_atomtype_nbparam(j,nf,atype);
+ c = sqrt(ci * cj);
plist->param[k].c[nf] = c;
}
+ }
+ }
break;
case eCOMB_ARITHMETIC:
/* c0 and c1 are epsilon and sigma */
for (i=k=0; (i < nr); i++)
for (j=0; (j < nr); j++,k++) {
- plist->param[k].c[0] = (atype->nb[i].C0+atype->nb[j].C0)*0.5;
- plist->param[k].c[1] = sqrt(atype->nb[i].C1*atype->nb[j].C1);
+ ci0 = get_atomtype_nbparam(i,0,atype);
+ cj0 = get_atomtype_nbparam(j,0,atype);
+ ci1 = get_atomtype_nbparam(i,1,atype);
+ cj1 = get_atomtype_nbparam(j,1,atype);
+ plist->param[k].c[0] = (ci0+cj0)*0.5;
+ plist->param[k].c[1] = sqrt(ci1*cj1);
}
break;
/* c0 and c1 are epsilon and sigma */
for (i=k=0; (i < nr); i++)
for (j=0; (j < nr); j++,k++) {
- plist->param[k].c[0] = sqrt(atype->nb[i].C0*atype->nb[j].C0);
- plist->param[k].c[1] = sqrt(atype->nb[i].C1*atype->nb[j].C1);
+ ci0 = get_atomtype_nbparam(i,0,atype);
+ cj0 = get_atomtype_nbparam(j,0,atype);
+ ci1 = get_atomtype_nbparam(i,1,atype);
+ cj1 = get_atomtype_nbparam(j,1,atype);
+ plist->param[k].c[0] = sqrt(ci0*cj0);
+ plist->param[k].c[1] = sqrt(ci1*ci1);
}
break;
/* Buckingham rules */
for(i=k=0; (i<nr); i++)
for(j=0; (j<nr); j++,k++) {
- plist->param[k].c[0] = sqrt(atype->nb[i].c[0] * atype->nb[j].c[0]);
- bi = atype->nb[i].c[1];
- bj = atype->nb[j].c[1];
+ ci0 = get_atomtype_nbparam(i,0,atype);
+ cj0 = get_atomtype_nbparam(j,0,atype);
+ ci2 = get_atomtype_nbparam(i,2,atype);
+ cj2 = get_atomtype_nbparam(j,2,atype);
+ bi = get_atomtype_nbparam(i,1,atype);
+ bj = get_atomtype_nbparam(j,1,atype);
+ plist->param[k].c[0] = sqrt(ci0 * cj0);
if ((bi == 0) || (bj == 0))
plist->param[k].c[1] = 0;
else
plist->param[k].c[1] = 2.0/(1/bi+1/bj);
- plist->param[k].c[2] = sqrt(atype->nb[i].c[2] * atype->nb[j].c[2]);
+ plist->param[k].c[2] = sqrt(ci2 * cj2);
}
break;
}
}
-void push_at (t_symtab *symtab, t_atomtype *at, t_bond_atomtype *bat,
+void push_at (t_symtab *symtab, t_atomtype at, t_bond_atomtype bat,
char *line,int nb_funct,
t_nbparam ***nbparam,t_nbparam ***pair)
{
double radius,vol,surftens;
char tmpfield[12][100]; /* Max 12 fields of width 100 */
char errbuf[256];
+ t_atom *atom;
+ t_param *param;
int atomnr;
bool have_atomic_number;
bool have_bonded_type;
+ snew(atom,1);
+ snew(param,1);
+
/* First assign input line to temporary array */
nfields=sscanf(line,"%s%s%s%s%s%s%s%s%s%s%s%s",
tmpfield[0],tmpfield[1],tmpfield[2],tmpfield[3],tmpfield[4],tmpfield[5],
pt=xl[j].ptype;
if (debug)
fprintf(debug,"ptype: %s\n",ptype_str[pt]);
-
- nr = bat->nr;
- /* First test if we have a new bond_atomtype */
- for (i=0; ((i<nr) && (strcasecmp(*(bat->atomname[i]),btype) != 0)); i++)
- ;
- if ((i==nr) || (nr==0)) {
- /* New bond_atomtype */
- srenew(bat->atomname,nr+1);
- bat->nr++;
- bat->atomname[nr] = put_symtab(symtab,btype);
- }
- batype_nr=i;
- nr = at->nr;
- /* Test if this atomtype overwrites another */
- for (i=0; ((i<nr) && (strcasecmp(*(at->atomname[i]),type) != 0)); i++)
- ;
-
- if ((i==nr) || (nr==0)) {
- /* New atomtype, get new space for arrays */
- srenew(at->atom,nr+1);
- srenew(at->atomname,nr+1);
- srenew(at->bondatomtype,nr+1);
- srenew(at->nb,nr+1);
- srenew(at->radius,nr+1);
- srenew(at->vol,nr+1);
- srenew(at->surftens,nr+1);
- srenew(at->atomnumber,nr+1);
- at->nr++;
+ atom->q = q;
+ atom->m = m;
+ atom->ptype = pt;
+ for (i=0; (i<MAXFORCEPARAM); i++)
+ param->c[i] = c[i];
+
+ if ((batype_nr = get_bond_atomtype_type(btype,bat)) == NOTSET)
+ add_bond_atomtype(bat,symtab,btype);
+ batype_nr = get_bond_atomtype_type(btype,bat);
+
+ if ((nr = get_atomtype_type(type,at)) != NOTSET) {
+ sprintf(errbuf,"Overriding atomtype %s",type);
+ warning(errbuf);
+ if ((nr = set_atomtype(nr,at,symtab,atom,type,param,batype_nr,
+ radius,vol,surftens,atomnr)) == NOTSET)
+ gmx_fatal(FARGS,"Replacing atomtype %s failed",type);
+ }
+ else if ((nr = add_atomtype(at,symtab,atom,type,param,
+ batype_nr,radius,vol,
+ surftens,atomnr)) == NOTSET)
+ gmx_fatal(FARGS,"Adding atomtype %s failed",type);
+ else {
/* Add space in the non-bonded parameters matrix */
- srenew(*nbparam,at->nr);
- snew((*nbparam)[at->nr-1],at->nr);
+ int atnr = get_atomtype_ntypes(at);
+ srenew(*nbparam,atnr);
+ snew((*nbparam)[atnr-1],atnr);
if (pair) {
- srenew(*pair,at->nr);
- snew((*pair)[at->nr-1],at->nr);
+ srenew(*pair,atnr);
+ snew((*pair)[atnr-1],atnr);
}
}
- else {
- sprintf(errbuf,"Overriding atomtype %s",type);
- warning(errbuf);
- nr = i;
- }
-
- /* fill the arrays */
- at->atomname[nr] = put_symtab(symtab,type);
- at->bondatomtype[nr] = batype_nr;
- at->atom[nr].ptype = pt;
- at->atom[nr].m = m;
- at->atom[nr].q = q;
- at->radius[nr] = radius;
- at->vol[nr] = vol;
- at->surftens[nr] = surftens;
- at->atomnumber[nr] = atomnr;
-
- for (i=0; (i<MAXFORCEPARAM); i++)
- at->nb[nr].c[i] = c[i];
+ sfree(atom);
+ sfree(param);
}
static void push_bondtype(t_params *bt,t_param *b,int nral,int ftype,
}
}
-void push_bt(directive d,t_params bt[],int nral,char ***typenames, int ntypes,char *line)
+void push_bt(directive d,t_params bt[],int nral,
+ t_atomtype at,
+ t_bond_atomtype bat,char *line)
{
const char *formal[MAXATOMLIST+1] = {
"%s",
t_param p;
char errbuf[256];
+ if ((bat && at) || (!bat && !at))
+ gmx_incons("You should pass either bat or at to push_bt");
+
/* Make format string (nral ints+functype) */
if ((nn=sscanf(line,formal[nral],
alc[0],alc[1],alc[2],alc[3],alc[4],alc[5])) != nral+1) {
c[i] = 0.0;
}
}
- for(i=0; (i<nral); i++)
- p.a[i]=name2index(alc[i],typenames,ntypes);
+ for(i=0; (i<nral); i++) {
+ if (at && ((p.a[i]=get_atomtype_type(alc[i],at)) == NOTSET))
+ gmx_fatal(FARGS,"Unknown atomtype %s\n",alc[i]);
+ else if (bat && ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET))
+ gmx_fatal(FARGS,"Unknown bond_atomtype %s\n",alc[i]);
+ }
for(i=0; (i<nrfp); i++)
p.c[i]=c[i];
push_bondtype (&(bt[ftype]),&p,nral,ftype,line);
}
-void push_dihedraltype(directive d,t_params bt[],char ***typenames,int ntypes,char *line)
+void push_dihedraltype(directive d,t_params bt[],
+ t_bond_atomtype bat,char *line)
{
const char *formal[MAXATOMLIST+1] = {
"%s",
for(i=0; (i<4); i++) {
if(!strcmp(alc[i],"X"))
p.a[i]=-1;
- else
- p.a[i]=name2index(alc[i],typenames,ntypes);
+ else {
+ if ((p.a[i]=get_bond_atomtype_type(alc[i],bat)) == NOTSET)
+ gmx_fatal(FARGS,"Unknown bond_atomtype %s",alc[i]);
+ }
}
for(i=0; (i<nrfp); i++)
p.c[i]=c[i];
}
-void push_nbt(directive d,t_nbparam **nbt,t_atomtype *atype,
+void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
char *pline,int nb_funct)
{
/* swap the atoms */
cr[i] = c[i];
/* Put the parameters in the matrix */
- ai = at2type (a0,atype);
- aj = at2type (a1,atype);
+ if ((ai = get_atomtype_type (a0,atype)) == NOTSET)
+ gmx_fatal(FARGS,"Atomtype %s not found",a0);
+ if ((aj = get_atomtype_type (a1,atype)) == NOTSET)
+ gmx_fatal(FARGS,"Atomtype %s not found",a1);
nbp = &(nbt[max(ai,aj)][min(ai,aj)]);
if (nbp->bSet) {
}
void push_atom(t_symtab *symtab,t_block *cgs,
- t_atoms *at,t_atomtype *atype,char *line,int *lastcg)
+ t_atoms *at,t_atomtype atype,char *line,int *lastcg)
{
int nr,ptype;
int resnumber,cgnumber,atomnr,type,typeB,nscan;
return;
}
sscanf(id,"%d",&atomnr);
- type = at2type(ctype,atype);
- ptype = atype->atom[type].ptype;
+ if ((type = get_atomtype_type(ctype,atype)) == NOTSET)
+ gmx_fatal(FARGS,"Atomtype %s not found",ctype);
+ ptype = get_atomtype_ptype(type,atype);
/* Set default from type */
- q0 = atype->atom[type].q;
- m0 = atype->atom[type].m;
+ q0 = get_atomtype_qA(type,atype);
+ m0 = get_atomtype_massA(type,atype);
typeB = type;
qB = q0;
mB = m0;
/* Optional parameters */
- nscan=sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
- &q,&m,ctypeB,&qb,&mb,check);
+ nscan = sscanf(line,"%*s%*s%*s%*s%*s%*s%lf%lf%s%lf%lf%s",
+ &q,&m,ctypeB,&qb,&mb,check);
/* Nasty switch that falls thru all the way down! */
if (nscan > 0) {
if (nscan > 1) {
m0 = mB = m;
if (nscan > 2) {
- typeB=at2type(ctypeB,atype);
- qB = atype->atom[typeB].q;
- mB = atype->atom[typeB].m;
+ typeB=get_atomtype_type(ctypeB,atype);
+ qB = get_atomtype_qA(typeB,atype);
+ mB = get_atomtype_massA(typeB,atype);
if (nscan > 3) {
qB = qb;
if (nscan > 4) {
static bool default_params(int ftype,t_params bt[],
- t_atoms *at,t_atomtype *atype,
+ t_atoms *at,t_atomtype atype,
t_param *p,bool bB,
t_param **param_def)
{
int nral = NRAL(ftype);
int nrfpA= interaction_function[ftype].nrfpA;
int nrfpB= interaction_function[ftype].nrfpB;
- int *batype=atype->bondatomtype;
if ((!bB && nrfpA == 0) || (bB && nrfpB == 0))
return TRUE;
if (bB) {
for (i=0; ((i < nr) && !bFound); i++) {
pi=&(bt[ftype].param[i]);
- bFound=( ((pi->AI==-1) || (batype[at->atom[p->AI].typeB]==pi->AI)) &&
- ((pi->AJ==-1) || (batype[at->atom[p->AJ].typeB]==pi->AJ)) &&
- ((pi->AK==-1) || (batype[at->atom[p->AK].typeB]==pi->AK)) &&
- ((pi->AL==-1) || (batype[at->atom[p->AL].typeB]==pi->AL)) );
+ bFound=
+ (((pi->AI==-1) ||
+ (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
+ ((pi->AJ==-1) ||
+ (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
+ ((pi->AK==-1) ||
+ (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
+ ((pi->AL==-1) ||
+ (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
+ );
}
} else {
for (i=0; ((i < nr) && !bFound); i++) {
pi=&(bt[ftype].param[i]);
- bFound=( ((pi->AI==-1) || (batype[at->atom[p->AI].type]==pi->AI)) &&
- ((pi->AJ==-1) || (batype[at->atom[p->AJ].type]==pi->AJ)) &&
- ((pi->AK==-1) || (batype[at->atom[p->AK].type]==pi->AK)) &&
- ((pi->AL==-1) || (batype[at->atom[p->AL].type]==pi->AL)) );
+ bFound=
+ (((pi->AI==-1) ||
+ (get_atomtype_batype(at->atom[p->AI].typeB,atype)==pi->AI)) &&
+ ((pi->AJ==-1) ||
+ (get_atomtype_batype(at->atom[p->AJ].typeB,atype)==pi->AJ)) &&
+ ((pi->AK==-1) ||
+ (get_atomtype_batype(at->atom[p->AK].typeB,atype)==pi->AK)) &&
+ ((pi->AL==-1) ||
+ (get_atomtype_batype(at->atom[p->AL].typeB,atype)==pi->AL))
+ );
}
}
} else { /* Not a dihedral */
pi=&(bt[ftype].param[i]);
if (bB)
for (j=0; ((j < nral) &&
- (atype->bondatomtype[at->atom[p->a[j]].typeB] == pi->a[j])); j++);
+ (get_atomtype_batype(at->atom[p->a[j]].typeB,atype)==pi->a[j])); j++)
+ ;
else
for (j=0; ((j < nral) &&
- (atype->bondatomtype[at->atom[p->a[j]].type] == pi->a[j])); j++);
+ (get_atomtype_batype(at->atom[p->a[j]].type,atype)==pi->a[j])); j++)
+ ;
bFound=(j==nral);
}
}
}
void push_bond(directive d,t_params bondtype[],t_params bond[],
- t_atoms *at,t_atomtype *atype,char *line,
+ t_atoms *at,t_atomtype atype,char *line,
bool bBonded,bool bGenPairs,bool bZero,bool *bWarn_copy_A_B)
{
const char *aaformat[MAXATOMLIST]= {
}
void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
- t_atoms *at,t_atomtype *atype,char *line)
+ t_atoms *at,t_atomtype atype,char *line)
{
char *ptr;
int type,ftype,j,n,ret,nj,a;
#include "typedefs.h"
#include "toputil.h"
+#include "gpp_atomtype.h"
+#include "gpp_bond_atomtype.h"
typedef struct {
int nr; /* The number of entries in the list */
} t_block2;
extern void generate_nbparams(int comb,int funct,t_params plist[],
- t_atomtype *atype);
+ t_atomtype atype);
-extern void push_at (t_symtab *symtab,t_atomtype *at,t_bond_atomtype *bat,char *line,int nb_funct,
+extern void push_at (t_symtab *symtab,t_atomtype at,
+ t_bond_atomtype bat,char *line,int nb_funct,
t_nbparam ***nbparam,t_nbparam ***pair);
-extern void push_bt(directive d,t_params bt[], int nral,
- char ***typenames, int ntypes, char *line);
+extern void push_bt(directive d,t_params bt[], int nral,
+ t_atomtype at,t_bond_atomtype bat,char *line);
extern void push_dihedraltype(directive d,t_params bt[],
- char ***typenames, int ntypes,char *line);
+ t_bond_atomtype bat,char *line);
-extern void push_nbt(directive d,t_nbparam **nbt,t_atomtype *atype,
+extern void push_nbt(directive d,t_nbparam **nbt,t_atomtype atype,
char *plines,int nb_funct);
extern void push_atom(t_symtab *symtab,
t_block *cgs,
t_atoms *at,
- t_atomtype *atype,
+ t_atomtype atype,
char *line,
int *lastcg);
extern void push_bondnow (t_params *bond, t_param *b);
extern void push_bond(directive d,t_params bondtype[],t_params bond[],
- t_atoms *at,t_atomtype *atype,char *line,
+ t_atoms *at,t_atomtype atype,char *line,
bool bBonded,bool bGenPairs,
bool bZero,bool *bWarn_copy_A_B);
extern void push_vsitesn(directive d,t_params bondtype[],t_params bond[],
- t_atoms *at,t_atomtype *atype,char *line);
+ t_atoms *at,t_atomtype atype,char *line);
extern void push_mol(int nrmols,t_molinfo mols[],char *pline,
int *whichmol,int *nrcopies);
return nh;
}
-void make_shake (t_params plist[],t_atoms *atoms,t_atomtype *at,int nshake)
+void make_shake (t_params plist[],t_atoms *atoms,t_atomtype at,int nshake)
{
char ***info=atoms->atomname;
t_params *pr;
#include "topio.h"
-void make_shake (t_params plist[],t_atoms *atoms,t_atomtype *at,int nshake);
+void make_shake (t_params plist[],t_atoms *atoms,t_atomtype at,int nshake);
#endif /* _topshake_h */
#include <config.h>
#endif
+#include <math.h>
#include "smalloc.h"
#include "sysstuff.h"
#include "macros.h"
#include "toputil.h"
#include "symtab.h"
#include "gmx_fatal.h"
-#include <math.h>
+#include "gpp_atomtype.h"
/* UTILITIES */
-int at2type(char *str,t_atomtype *at)
-{
- int i;
-
- for (i=0; (i<at->nr) && strcasecmp(str,*(at->atomname[i])); i++)
- ;
- if (i == at->nr)
- gmx_fatal(FARGS,"Atomtype '%s' not found!",str);
-
- return i;
-}
-
-
-int name2index(char *str, char ***typenames, int ntypes)
-{
- int i;
-
- for (i=0; (i<ntypes) && strcasecmp(str,*(typenames[i])); i++)
- ;
- if (i == ntypes)
- gmx_fatal(FARGS,"Bonded/nonbonded atom type '%s' not found!",str);
-
- return i;
-}
-
-
-char *type2nm(int nt, t_atomtype *at)
-{
- if ((nt < 0) || (nt >= at->nr))
- gmx_fatal(FARGS,"nt out of range: %d",nt);
-
- return *(at->atomname[nt]);
-}
-
void set_p_string(t_param *p,char *s)
{
if (s) {
}
}
-/* INIT STRUCTURES */
-
-void init_atomtype (t_atomtype *at)
-{
- at->nr = 0;
- at->atom = NULL;
- at->atomname = NULL;
- at->nb = NULL;
- at->bondatomtype = NULL;
- at->radius = NULL;
- at->vol = NULL;
- at->surftens = NULL;
- at->atomnumber = NULL;
-}
-
-void init_bond_atomtype (t_bond_atomtype *bat)
-{
- bat->nr = 0;
- bat->atomname = NULL;
-}
-
-
void init_plist(t_params plist[])
{
int i;
/* PRINTING STRUCTURES */
-void print_at (FILE * out, t_atomtype *at)
-{
- int i;
- t_atom *atom = at->atom;
- t_param *nb = at->nb;
-
- fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
- fprintf (out,"; %6s %8s %8s %8s %12s %12s\n",
- "type","mass","charge","particle","c6","c12");
- for (i=0; (i<at->nr); i++)
- fprintf(out,"%8s %8.3f %8.3f %8s %12e %12e\n",
- *(at->atomname[i]),atom[i].m,atom[i].q,"A",
- nb[i].C0,nb[i].C1);
-
- fprintf (out,"\n");
-}
-
-static void print_nbt (FILE *out,char *title,t_atomtype *at,
+static void print_nbt (FILE *out,char *title,t_atomtype at,
int ftype,t_params *nbt)
{
- int f,i,j,k,l,nrfp;
+ int f,i,j,k,l,nrfp,ntype;
if (ftype == F_LJ)
f=1;
fprintf (out,"\n");
/* print non-bondtypes */
- for (i=k=0; (i<at->nr); i++)
- for(j=0; (j<at->nr); j++,k++) {
+ ntype = get_atomtype_ntypes(at);
+ for (i=k=0; (i<ntype); i++)
+ for(j=0; (j<ntype); j++,k++) {
fprintf (out,"%8s %8s %8d",
- *(at->atomname[i]),*(at->atomname[j]),f);
+ get_atomtype_name(i,at),get_atomtype_name(f,at),f);
for(l=0; (l<nrfp); l++)
fprintf (out," %12.5e",nbt->param[k].c[l]);
fprintf (out,"\n");
fprintf (out,"\n");
}
-void print_bt(FILE *out, directive d, t_atomtype *at,
+void print_bt(FILE *out, directive d, t_atomtype at,
int ftype,int fsubtype,t_params plist[],
bool bFullDih)
{
bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
if (!bDih)
for (j=0; (j<nral); j++)
- fprintf (out,"%5s ",*(at->atomname[bt->param[i].a[j]]));
+ fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[j],at));
else
for(j=0; (j<2); j++)
- fprintf (out,"%5s ",*(at->atomname[bt->param[i].a[dihp[f][j]]]));
+ fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[dihp[f][j]],at));
fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
if (bt->param[i].s[0])
}
}
-void print_atoms(FILE *out,t_atomtype *atype,t_atoms *at,int *cgnr)
+void print_atoms(FILE *out,t_atomtype atype,t_atoms *at,int *cgnr)
{
int i;
- int itype;
- char *as;
+ int tpA,tpB;
+ char *as,*tpnmA,*tpnmB;
double qtot;
as=dir2str(d_atoms);
if (at->nres) {
/* if the information is present... */
for (i=0; (i < at->nr); i++) {
- itype=at->atom[i].type;
- if ((itype < 0) || (itype > atype->nr))
- gmx_fatal(FARGS,"itype = %d, i= %d in print_atoms",itype,i);
-
+ tpA = at->atom[i].type;
+ if ((tpnmA = get_atomtype_name(tpA,atype)) == NULL)
+ gmx_fatal(FARGS,"tpA = %d, i= %d in print_atoms",tpA,i);
+
fprintf(out,"%6d %10s %6d %6s %6s %6d %10g %10g",
- i+1,*(atype->atomname[itype]),
- at->atom[i].resnr+1,
+ i+1,tpnmA,at->atom[i].resnr+1,
*(at->resname[at->atom[i].resnr]),
*(at->atomname[i]),cgnr[i],
at->atom[i].q,at->atom[i].m);
if (PERTURBED(at->atom[i])) {
+ tpB = at->atom[i].typeB;
+ if ((tpnmB = get_atomtype_name(tpB,atype)) == NULL)
+ gmx_fatal(FARGS,"tpB = %d, i= %d in print_atoms",tpB,i);
fprintf(out," %6s %10g %10g",
- *(atype->atomname[at->atom[i].typeB]),
- at->atom[i].qB,at->atom[i].mB);
+ tpnmB,at->atom[i].qB,at->atom[i].mB);
}
qtot += (double)at->atom[i].q;
if ( fabs(qtot) < 4*GMX_REAL_EPS )
{
t_symtab stab;
t_atomtype atype;
+ t_param *param;
+ t_atom *a;
int i;
- snew(atype.atom,natoms);
- snew(atype.atomname,natoms);
+ atype = init_atomtype();
+ snew(a,1);
+ snew(param,1);
open_symtab(&stab);
for (i=0; (i < natoms); i++) {
char buf[12];
sprintf(buf,"%4d",(i+1));
- atype.atomname[i]=put_symtab(&stab,buf);
+ add_atomtype(atype,&stab,a,buf,param,0,0,0,0,0);
}
- print_bt(out,d,&atype,ftype,fsubtype,plist,TRUE);
+ print_bt(out,d,atype,ftype,fsubtype,plist,TRUE);
done_symtab(&stab);
- sfree(atype.atom);
- sfree(atype.atomname);
+ sfree(a);
+ sfree(param);
+ done_atomtype(&atype);
}
#define _toputil_h
#include "grompp.h"
+#include "gpp_atomtype.h"
/* UTILITIES */
-extern int at2type(char *str, t_atomtype *at);
-
extern int name2index(char *str, char ***typenames, int ntypes);
-extern char *type2nm(int nt, t_atomtype *at);
-
extern void pr_alloc (int extra, t_params *pr);
extern void set_p_string(t_param *p,char *s);
extern void init_plist(t_params plist[]);
-extern void init_atomtype (t_atomtype *at);
-
-extern void init_bond_atomtype (t_bond_atomtype *bat);
-
extern void init_molinfo(t_molinfo *mol);
extern void init_top (t_topology *top);
extern void print_blocka(FILE *out,char *szName,char *szIndex,
char *szA,t_blocka *block);
-extern void print_atoms(FILE *out,t_atomtype *atype,t_atoms *at,int *cgnr);
+extern void print_atoms(FILE *out,t_atomtype atype,t_atoms *at,int *cgnr);
extern void print_bondeds(FILE *out,int natoms,directive d,
int ftype,int fsubtype,t_params plist[]);
return angle;
}
-static bool calc_vsite3_param(t_atomtype *atype,
+static bool calc_vsite3_param(t_atomtype atype,
t_param *param, t_atoms *at,
int nrbond, t_mybonded *bonds,
int nrang, t_mybonded *angles )
int i;
for (i=0; i<4; i++)
fprintf(debug,"atom %u type %s ",
- param->a[i]+1,type2nm(at->atom[param->a[i]].type,atype));
+ param->a[i]+1,
+ get_atomtype_name(at->atom[param->a[i]].type,atype));
fprintf(debug,"\n");
}
bXH3 =
- ( (strncasecmp(type2nm(at->atom[param->AK].type,atype),"MNH",3)==0) &&
- (strncasecmp(type2nm(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
- ( (strncasecmp(type2nm(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
- (strncasecmp(type2nm(at->atom[param->AL].type,atype),"MCH3",4)==0) );
+ ( (strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MNH",3)==0) &&
+ (strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
+ ( (strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
+ (strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MCH3",4)==0) );
bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
return bError;
}
-static bool calc_vsite3out_param(t_atomtype *atype,
+static bool calc_vsite3out_param(t_atomtype atype,
t_param *param, t_atoms *at,
int nrbond, t_mybonded *bonds,
int nrang, t_mybonded *angles)
int i;
for (i=0; i<4; i++)
fprintf(debug,"atom %u type %s ",
- param->a[i]+1,type2nm(at->atom[param->a[i]].type,atype));
+ param->a[i]+1,get_atomtype_name(at->atom[param->a[i]].type,atype));
fprintf(debug,"\n");
}
bXH3 =
- ( (strncasecmp(type2nm(at->atom[param->AK].type,atype),"MNH",3)==0) &&
- (strncasecmp(type2nm(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
- ( (strncasecmp(type2nm(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
- (strncasecmp(type2nm(at->atom[param->AL].type,atype),"MCH3",4)==0) );
+ ( (strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MNH",3)==0) &&
+ (strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MNH",3)==0) ) ||
+ ( (strncasecmp(get_atomtype_name(at->atom[param->AK].type,atype),"MCH3",4)==0) &&
+ (strncasecmp(get_atomtype_name(at->atom[param->AL].type,atype),"MCH3",4)==0) );
/* check if construction parity must be swapped */
bSwapParity = ( param->C1 == -1 );
-int set_vsites(bool bVerbose, t_atoms *atoms, t_atomtype *atype,
+int set_vsites(bool bVerbose, t_atoms *atoms, t_atomtype atype,
t_params plist[])
{
int i,j,ftype;
#include "typedefs.h"
#include "grompp.h"
+#include "gpp_atomtype.h"
-extern int set_vsites(bool bVerbose, t_atoms *atoms, t_atomtype *atype,
+extern int set_vsites(bool bVerbose, t_atoms *atoms, t_atomtype atype,
t_params plist[]);
/* set parameters for vritual sites, return number of virtual sites */
t_params plist[F_NRE];
t_excls *excls;
t_atoms *atoms; /* list with all atoms */
- t_atomtype *atype;
+ t_atomtype atype;
t_nextnb nnb;
x2top_nm2t nm2t;
x2top_qat qa;
static char *molnm = "ICE";
static char *ff = "select";
static char *qgen[] = { NULL, "None", "Linear", "Yang", "Bultinck",
- "SM1", "SM2", "SM3", "SM4", NULL };
+ "SMp", "SMs", "SMps", "SMg", "SMgs", NULL };
t_pargs pa[] = {
{ "-ff", FALSE, etSTR, {&ff},
"Select the force field for your simulation." },
mk_bonds(nm2t,atoms,x,&(plist[F_BONDS]),nbonds,forcefield,
bPBC,box,atomprop,btol);
+ /* Setting the atom types: this depends on the bonding */
open_symtab(&symtab);
atype = set_atom_type(&symtab,atoms,&(plist[F_BONDS]),nbonds,nm2t);
- /* Read charges */
+ /* Read charges and polarizabilities, if provided */
if (opt2bSet("-d",NFILE,fnm))
qa = rd_q_alpha(opt2fn("-d",NFILE,fnm));
+ /* Check which algorithm to use for charge generation */
alg = eqgNone;
if (qgen[0]) {
alg = name2eemtype(qgen[0]);
maxiter,atomprop,qtotref);
if (bPolarize)
- add_shells(nm2t,atoms,atype,&(plist[F_BONDS]),&(plist[F_POLARIZATION]),
- &x,&symtab);
-
+ add_shells(nm2t,&atoms,atype,&(plist[F_BONDS]),
+ &(plist[F_POLARIZATION]),&x,&symtab);
+
/* Make Angles and Dihedrals */
snew(excls,atoms->nr);
printf("Generating angles and dihedrals from bonds...\n");
print_nnb(&nnb,"NNB");
gen_pad(&nnb,atoms,bH14,nexcl,plist,excls,NULL,
bAllDih,bRemoveDih,TRUE);
- delete_shell_interactions(plist,atoms,atype,&nnb,excls);
+ /*delete_shell_interactions(plist,atoms,atype,&nnb,excls);*/
done_nnb(&nnb);
mu = calc_dip(atoms,x);
#include "grompp.h"
#include "add_par.h"
#include "gmx_random.h"
+#include "gpp_atomtype.h"
void mk_bonds(x2top_nm2t nmt,
t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
return cgnr;
}
-t_atomtype *set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
- int *nbonds,x2top_nm2t nm2t)
+t_atomtype set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
+ int *nbonds,x2top_nm2t nm2t)
{
- t_atomtype *atype;
+ t_atomtype atype;
int nresolved;
- snew(atype,1);
+ atype = init_atomtype();
snew(atoms->atomtype,atoms->nr);
nresolved = nm2type(nm2t,tab,atoms,atype,nbonds,bonds);
if (nresolved != atoms->nr)
nresolved,atoms->nr);
fprintf(stderr,"There are %d different atom types in your sample\n",
- atype->nr);
+ get_atomtype_ntypes(atype));
return atype;
}
}
void delete_shell_interactions(t_params plist[F_NRE],t_atoms *atoms,
- t_atomtype *atype,t_nextnb *nnb,
+ t_atomtype atype,t_nextnb *nnb,
t_excls excls[])
{
int atp,jtp,jid,i,j,k,l,m,ftype,nb,nra,npol=0;
for(j=0; (j<plist[ftype].nr); j++) {
for(k=0; (k<nra); k++) {
atp = atoms->atom[p[j].a[k]].type;
- if (strcasecmp(*atype->atomname[atp],"SHELL") == 0) {
+ if (get_atomtype_ptype(atp,atype) == eptShell) {
bRemove[j] = TRUE;
if (ftype == F_BONDS) {
memcpy(&plist[F_POLARIZATION].param[npol],
/* now for all atoms */
for (i=0; (i < atoms->nr); i++) {
atp = atoms->atom[i].type;
- if (strcasecmp(*atype->atomname[atp],"SHELL") == 0) {
-
+ if (get_atomtype_ptype(atp,atype) == eptShell) {
for(m=3; (m<=4); m++) {
/* for all fifth bonded atoms of atom i */
for (j=0; (j < nnb->nrexcl[i][m]); j++) {
/* store the 1st neighbour in nb */
nb = nnb->a[i][m][j];
jtp = atoms->atom[nb].type;
- if ((i != nb) && (strcasecmp(*atype->atomname[jtp],"SHELL") == 0)) {
+ if ((i != nb) && (strcasecmp(get_atomtype_name(jtp,atype),"SHELL") == 0)) {
srenew(excls[i].e,excls[i].nr+1);
excls[i].e[excls[i].nr++] = nb;
fprintf(stderr,"Excluding %d from %d\n",nb+1,i+1);
jid = nnb->a[i][1][j];
atp = atoms->atom[jid].type;
- if (strcasecmp(*atype->atomname[atp],"SHELL") == 0) {
+ if (get_atomtype_ptype(atp,atype) == eptShell) {
bHaveShell[i] = TRUE;
shell_index[i] = jid;
}
atoms->atom[i].qB = atoms->atom[i].q;
}
-void add_shells(x2top_nm2t nm2t,t_atoms *atoms,
- t_atomtype *atype,
+void add_shells(x2top_nm2t nm2t,t_atoms **atoms,
+ t_atomtype atype,
t_params *bonds,t_params *pols,
rvec **x,t_symtab *symtab)
{
int i,j,iat,shell,atp,ns=0;
int *renum;
- char **ptr;
char buf[32];
t_param p;
+ t_atom *shell_atom;
+ t_atoms *newa;
+ rvec *newx;
+ snew(shell_atom,1);
memset(&p,0,sizeof(p));
- snew(renum,atoms->nr*2);
- for(i=0; (i<atoms->nr); i++) {
- atp = atoms->atom[i].type;
+ snew(renum,(*atoms)->nr*2);
+ for(i=0; (i<(*atoms)->nr); i++) {
+ atp = (*atoms)->atom[i].type;
renum[i] = i+ns;
- if (atype->atom[atp].qB != 0) {
+ if (get_atomtype_qB(atp,atype) != 0) {
ns++;
p.AI = renum[i];
p.AJ = renum[i]+1;
- p.C0 = atype->atom[atp].qB;
+ p.C0 = get_atomtype_qB(atp,atype);
push_bondnow(pols,&p);
}
}
- shell = atype->nr - 1;
+ shell_atom->ptype = eptShell;
+ shell = add_atomtype(atype,symtab,shell_atom,"SHELL",&p,
+ 0,0,0,0,0);
+
if (ns > 0) {
- srenew(atoms->atom,atoms->nr+ns);
- srenew(atoms->atomname,atoms->nr+ns);
- srenew((*x),atoms->nr+ns);
- for(i=atoms->nr; (i<atoms->nr+ns); i++) {
- memset(&(atoms->atom[i]),0,sizeof(atoms->atom[i]));
- atoms->atomname[i] = NULL;
- clear_rvec((*x)[i]);
+ snew(newa,1);
+ init_t_atoms(newa,(*atoms)->nr+ns,TRUE);
+ newa->nres = (*atoms)->nres;
+ snew(newx,(*atoms)->nr+ns);
+
+ for(i=0; (i<(*atoms)->nr); i++) {
+ newa->atom[renum[i]] = (*atoms)->atom[i];
+ newa->atomname[renum[i]] = put_symtab(symtab,*(*atoms)->atomname[i]);
+ copy_rvec(newx[renum[i]],(*x)[i]);
}
- for(i=atoms->nr-1; (i>=0); i--) {
- atoms->atom[renum[i]] = atoms->atom[i];
- atoms->atomname[renum[i]] = atoms->atomname[i];
- copy_rvec((*x)[renum[i]],(*x)[i]);
+ for(i=0; (i<(*atoms)->nres); i++) {
+ newa->resname[i] = put_symtab(symtab,*(*atoms)->resname[i]);
}
- for(i=0; (i<atoms->nr); i++) {
+
+ for(i=0; (i<(*atoms)->nr); i++) {
iat = renum[i];
for(j=iat+1; (j<renum[i+1]); j++) {
- atoms->atom[j] = atoms->atom[iat];
- atoms->atom[iat].q = 0;
- atoms->atom[iat].qB = 0;
- atoms->atom[j].m = 0;
- atoms->atom[j].mB = 0;
- atoms->atom[j].type = shell;
- sprintf(buf,"S%s",*(atoms->atomname[iat]));
- snew(ptr,1);
- *ptr = strdup(buf);
- atoms->atomname = *ptr;/*put_symtab(symtab,ptr);*/
+ newa->atom[j] = (*atoms)->atom[iat];
+ newa->atom[iat].q = 0;
+ newa->atom[iat].qB = 0;
+ newa->atom[j].m = 0;
+ newa->atom[j].mB = 0;
+ newa->atom[j].type = shell;
+ newa->atom[j].resnr = (*atoms)->atom[i].resnr;
+ sprintf(buf,"S%s",*((*atoms)->atomname[i]));
+ newa->atomname[j] = put_symtab(symtab,buf);
copy_rvec((*x)[j],(*x)[iat]);
}
}
+ *atoms = newa;
+ *x = newx;
}
- atoms->nr += ns;
sfree(renum);
+ sfree(shell_atom);
}
extern real calc_dip(t_atoms *atoms,rvec x[]);
extern void delete_shell_interactions(t_params plist[F_NRE],t_atoms *atoms,
- t_atomtype *atype,t_nextnb *nnb,
+ t_atomtype atype,t_nextnb *nnb,
t_excls excls[]);
extern void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[]);
t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
bool bPBC,matrix box,void *atomprop,real tol);
-extern t_atomtype *set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
- int *nbonds,x2top_nm2t nm2t);
+extern t_atomtype set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
+ int *nbonds,x2top_nm2t nm2t);
-extern void add_shells(x2top_nm2t nm2t,t_atoms *atoms,
- t_atomtype *atype,
+extern void add_shells(x2top_nm2t nm2t,t_atoms **atoms,
+ t_atomtype atype,
t_params *bonds,t_params *pols,
rvec **x,t_symtab *symtab);
} t_eemrecord;
static char *eemtype_name[eqgNR] = {
- "None", "Linear", "Yang", "Bultinck", "SM1", "SM2", "SM3", "SM4"
+ "None", "Linear", "Yang", "Bultinck", "SMp", "SMs", "SMps", "SMg", "SMpg",
};
int name2eemtype(char *name)
else
*wj = (3/(4*er->eep[index].w));
}
- else if ((er->eep[index].eemtype == eqgSM2) ||
- (er->eep[index].eemtype == eqgSM3) ||
- (er->eep[index].eemtype == eqgSM4)) {
+ else if ((er->eep[index].eemtype == eqgSMs) ||
+ (er->eep[index].eemtype == eqgSMps) ||
+ (er->eep[index].eemtype == eqgSMg) ||
+ (er->eep[index].eemtype == eqgSMpg)) {
Z = er->eep[index].elem;
*wj = 1.0/er->eep[index].w;
/**wj = (Z-q)/(Z*er->eep[index].w);*/
void matrix_invert(FILE *fp,int n,double **a)
{
int i,j,m,lda,*ipiv,lwork,info;
- double **test,**id,*work;
+ double **test=NULL,**id,*work;
if (fp) {
fprintf(fp,"Inverting %d square matrix\n",n);
/* Dump the database for debugging. Can be reread by the program */
extern int nm2type(x2top_nm2t nm2t,t_symtab *tab,t_atoms *atoms,
- t_atomtype *atype,int *nbonds,t_params *bond);
+ t_atomtype atype,int *nbonds,t_params *bond);
/* Try to determine the atomtype (force field dependent) for the atoms
* with help of the bond list
*/
switch (eemtype) {
case eqgBultinck:
- case eqgSM1:
+ case eqgSMp:
eNN = coul_nucl_nucl(0,r);
break;
- case eqgSM2:
- case eqgSM3:
- case eqgSM4:
+ case eqgSMs:
if ((wi > 0) && (wj > 0)) {
wij = 2*(wi + wj)/(wi*wj);
eSS = coul_slater_slater(wij,r);
eNN = coul_nucl_nucl(0,r);
}
break;
+ case eqgSMps:
+ case eqgSMg:
+ case eqgSMpg:
case eqgYang:
default:
gmx_fatal(FARGS,"Can not treat algorithm %d yet in calc_jab",eemtype);
if (i != j) {
wj = qgen->wj[j];
qgen->Jab[i][j] = calc_jab(qgen->x[i],qgen->x[j],wi,wj,qgen->eemtype);
- if ((eemtype == eqgSM3) || (eemtype == eqgSM4))
- qgen->rhs[i] -= qgen->elem[j]*calc_j1(qgen->x[i],qgen->x[j],wi,wj,qgen->eemtype);
+ if ((eemtype == eqgSMps) || (eemtype == eqgSMpg))
+ qgen->rhs[i] -= qgen->elem[j]*calc_j1(qgen->x[i],qgen->x[j],
+ wi,wj,qgen->eemtype);
}
}
}
please_cite(stdout,"Bultinck2002a");
generate_charges_bultinck(eem,atoms,x,tol,maxiter,atomprop);
break;
- case eqgSM1:
- case eqgSM2:
- case eqgSM3:
- case eqgSM4:
+ case eqgSMp:
+ case eqgSMs:
+ case eqgSMps:
+ case eqgSMg:
+ case eqgSMpg:
(void) generate_charges_sm(debug,molname,eem,atoms,x,tol,maxiter,atomprop,
qtotref,eemtype);
break;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff