{
maxsize = max(maxsize,cgs->index[cg+1]-cgs->index[cg]);
}
-
- if (maxsize > 10)
+
+ if (maxsize > MAX_CHARGEGROUP_SIZE)
+ {
+ gmx_fatal(FARGS,"The largest charge group contains %d atoms. The maximum is %d.",maxsize,MAX_CHARGEGROUP_SIZE);
+ }
+ else if (maxsize > 10)
{
set_warning_line(wi,topfn,-1);
sprintf(warn_buf,
}
}
+static void check_bonds_timestep(gmx_mtop_t *mtop,double dt,warninp_t wi)
+{
+ /* This check is not intended to ensure accurate integration,
+ * rather it is to signal mistakes in the mdp settings.
+ * A common mistake is to forget to turn on constraints
+ * for MD after energy minimization with flexible bonds.
+ * This check can also detect too large time steps for flexible water
+ * models, but such errors will often be masked by the constraints
+ * mdp options, which turns flexible water into water with bond constraints,
+ * but without an angle constraint. Unfortunately such incorrect use
+ * of water models can not easily be detected without checking
+ * for specific model names.
+ *
+ * The stability limit of leap-frog or velocity verlet is 4.44 steps
+ * per oscillational period.
+ * But accurate bonds distributions are lost far before that limit.
+ * To allow relatively common schemes (although not common with Gromacs)
+ * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
+ * we set the note limit to 10.
+ */
+ int min_steps_warn=5;
+ int min_steps_note=10;
+ t_iparams *ip;
+ int molt;
+ gmx_moltype_t *moltype,*w_moltype;
+ t_atom *atom;
+ t_ilist *ilist,*ilb,*ilc,*ils;
+ int ftype;
+ int i,a1,a2,w_a1,w_a2,j;
+ real twopi2,limit2,fc,re,m1,m2,period2,w_period2;
+ gmx_bool bFound,bWater,bWarn;
+ char warn_buf[STRLEN];
+
+ ip = mtop->ffparams.iparams;
+
+ twopi2 = sqr(2*M_PI);
+
+ limit2 = sqr(min_steps_note*dt);
+
+ w_a1 = w_a2 = -1;
+ w_period2 = -1.0;
+
+ w_moltype = NULL;
+ for(molt=0; molt<mtop->nmoltype; molt++)
+ {
+ moltype = &mtop->moltype[molt];
+ atom = moltype->atoms.atom;
+ ilist = moltype->ilist;
+ ilc = &ilist[F_CONSTR];
+ ils = &ilist[F_SETTLE];
+ for(ftype=0; ftype<F_NRE; ftype++)
+ {
+ if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
+ {
+ continue;
+ }
+
+ ilb = &ilist[ftype];
+ for(i=0; i<ilb->nr; i+=3)
+ {
+ fc = ip[ilb->iatoms[i]].harmonic.krA;
+ re = ip[ilb->iatoms[i]].harmonic.rA;
+ if (ftype == F_G96BONDS)
+ {
+ /* Convert squared sqaure fc to harmonic fc */
+ fc = 2*fc*re;
+ }
+ a1 = ilb->iatoms[i+1];
+ a2 = ilb->iatoms[i+2];
+ m1 = atom[a1].m;
+ m2 = atom[a2].m;
+ if (fc > 0 && m1 > 0 && m2 > 0)
+ {
+ period2 = twopi2*m1*m2/((m1 + m2)*fc);
+ }
+ else
+ {
+ period2 = GMX_FLOAT_MAX;
+ }
+ if (debug)
+ {
+ fprintf(debug,"fc %g m1 %g m2 %g period %g\n",
+ fc,m1,m2,sqrt(period2));
+ }
+ if (period2 < limit2)
+ {
+ bFound = FALSE;
+ for(j=0; j<ilc->nr; j+=3)
+ {
+ if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) ||
+ (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1))
+ {
+ bFound = TRUE;
+ }
+ }
+ for(j=0; j<ils->nr; j+=2)
+ {
+ if ((a1 >= ils->iatoms[j+1] && a1 < ils->iatoms[j+1]+3) &&
+ (a2 >= ils->iatoms[j+1] && a2 < ils->iatoms[j+1]+3))
+ {
+ bFound = TRUE;
+ }
+ }
+ if (!bFound &&
+ (w_moltype == NULL || period2 < w_period2))
+ {
+ w_moltype = moltype;
+ w_a1 = a1;
+ w_a2 = a2;
+ w_period2 = period2;
+ }
+ }
+ }
+ }
+ }
+
+ if (w_moltype != NULL)
+ {
+ bWarn = (w_period2 < sqr(min_steps_warn*dt));
+ /* A check that would recognize most water models */
+ bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' &&
+ w_moltype->atoms.nr <= 5);
+ sprintf(warn_buf,"The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n"
+ "%s",
+ *w_moltype->name,
+ w_a1+1,*w_moltype->atoms.atomname[w_a1],
+ w_a2+1,*w_moltype->atoms.atomname[w_a2],
+ sqrt(w_period2),bWarn ? min_steps_warn : min_steps_note,dt,
+ bWater ?
+ "Maybe you asked for fexible water." :
+ "Maybe you forgot to change the constraints mdp option.");
+ if (bWarn)
+ {
+ warning(wi,warn_buf);
+ }
+ else
+ {
+ warning_note(wi,warn_buf);
+ }
+ }
+}
+
static void check_vel(gmx_mtop_t *mtop,rvec v[])
{
gmx_mtop_atomloop_all_t aloop;
return count;
}
+static void check_gbsa_params_charged(gmx_mtop_t *sys, gpp_atomtype_t atype)
+{
+ int i,nmiss,natoms,mt;
+ real q;
+ const t_atoms *atoms;
+
+ nmiss = 0;
+ for(mt=0;mt<sys->nmoltype;mt++)
+ {
+ atoms = &sys->moltype[mt].atoms;
+ natoms = atoms->nr;
+
+ for(i=0;i<natoms;i++)
+ {
+ q = atoms->atom[i].q;
+ if ((get_atomtype_radius(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_vol(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_surftens(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_gb_radius(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_S_hct(atoms->atom[i].type,atype) == 0) &&
+ q != 0)
+ {
+ fprintf(stderr,"\nGB parameter(s) zero for atom type '%s' while charge is %g\n",
+ get_atomtype_name(atoms->atom[i].type,atype),q);
+ nmiss++;
+ }
+ }
+ }
+
+ if (nmiss > 0)
+ {
+ gmx_fatal(FARGS,"Can't do GB electrostatics; the implicit_genborn_params section of the forcefield has parameters with value zero for %d atomtypes that occur as charged atoms.",nmiss);
+ }
+}
+
+
static void check_gbsa_params(t_inputrec *ir,gpp_atomtype_t atype)
{
int nmiss,i;
get_atomtype_gb_radius(i,atype) < 0 ||
get_atomtype_S_hct(i,atype) < 0)
{
- fprintf(stderr,"GB parameter(s) missing or negative for atom type '%s'\n",
+ fprintf(stderr,"\nGB parameter(s) missing or negative for atom type '%s'\n",
get_atomtype_name(i,atype));
nmiss++;
}
if (nmiss > 0)
{
- gmx_fatal(FARGS,"Can't do GB electrostatics; the forcefield is missing %d values for\n"
- "atomtype radii, or they might be negative\n.",nmiss);
+ gmx_fatal(FARGS,"Can't do GB electrostatics; the implicit_genborn_params section of the forcefield is missing parameters for %d atomtypes or they might be negative.",nmiss);
}
}
char warn_buf[STRLEN];
t_filenm fnm[] = {
- { efMDP, NULL, NULL, ffOPTRD },
+ { efMDP, NULL, NULL, ffREAD },
{ efMDP, "-po", "mdout", ffWRITE },
{ efSTX, "-c", NULL, ffREAD },
{ efSTX, "-r", NULL, ffOPTRD },
{
/* Now we have renumbered the atom types, we can check the GBSA params */
check_gbsa_params(ir,atype);
+
+ /* Check that all atoms that have charge and/or LJ-parameters also have
+ * sensible GB-parameters
+ */
+ check_gbsa_params_charged(sys,atype);
}
/* PELA: Copy the atomtype data to the topology atomtype list */
check_cg_sizes(ftp2fn(efTOP,NFILE,fnm),&sys->moltype[i].cgs,wi);
}
+ if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
+ {
+ check_bonds_timestep(sys,ir->delta_t,wi);
+ }
+
check_warning_error(wi,FARGS);
if (bVerbose)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff