wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
static char foreign_lambda[STRLEN];
static char **pull_grp;
+static char **rot_grp;
static char anneal[STRLEN],anneal_npoints[STRLEN],
anneal_time[STRLEN],anneal_temp[STRLEN];
static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
{
dt_pcoupl = ir->nstpcouple*ir->delta_t;
- sprintf(err_buf,"tau_p must be > 0 instead of %g\n",ir->tau_p);
+ sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
CHECK(ir->tau_p <= 0);
if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
{
- sprintf(warn_buf,"For proper integration of the %s barostat, tau_p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
+ sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
warning(wi,warn_buf);
}
sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
CHECK(ir->coulombtype == eelPPPM);
-
+
+ if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
+ {
+ sprintf(warn_buf,
+ "You are generating velocities so I am assuming you "
+ "are equilibrating a system. You are using "
+ "Parrinello-Rahman pressure coupling, but this can be "
+ "unstable for equilibration. If your system crashes, try "
+ "equilibrating first with Berendsen pressure coupling. If "
+ "you are not equilibrating the system, you can probably "
+ "ignore this warning.");
+ warning(wi,warn_buf);
+ }
}
else if (ir->coulombtype == eelPPPM)
{
}
if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
- sprintf(warn_buf,"epsilon_r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
+ sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
warning_note(wi,warn_buf);
}
if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
- sprintf(warn_buf,"epsilon_r = %g and epsilon_rf = 1 with reaction field, assuming old format and exchanging epsilon_r and epsilon_rf",ir->epsilon_r);
+ sprintf(warn_buf,"epsilon-r = %g and epsilon-rf = 1 with reaction field, assuming old format and exchanging epsilon-r and epsilon-rf",ir->epsilon_r);
warning(wi,warn_buf);
ir->epsilon_rf = ir->epsilon_r;
ir->epsilon_r = 1.0;
}
if (getenv("GALACTIC_DYNAMICS") == NULL) {
- sprintf(err_buf,"epsilon_r must be >= 0 instead of %g\n",ir->epsilon_r);
+ sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
CHECK(ir->epsilon_r < 0);
}
/* reaction field (at the cut-off) */
if (ir->coulombtype == eelRF_ZERO) {
- sprintf(err_buf,"With coulombtype = %s, epsilon_rf must be 0",
+ sprintf(err_buf,"With coulombtype = %s, epsilon-rf must be 0",
eel_names[ir->coulombtype]);
CHECK(ir->epsilon_rf != 0);
}
- sprintf(err_buf,"epsilon_rf must be >= epsilon_r");
+ sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
(ir->epsilon_r == 0));
if (ir->epsilon_rf == ir->epsilon_r) {
- sprintf(warn_buf,"Using epsilon_rf = epsilon_r with %s does not make sense",
+ sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
eel_names[ir->coulombtype]);
warning(wi,warn_buf);
}
if (EEL_PME(ir->coulombtype)) {
if (ir->pme_order < 3) {
- warning_error(wi,"pme_order can not be smaller than 3");
+ warning_error(wi,"pme-order can not be smaller than 3");
}
}
if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
if (ir->ewald_geometry == eewg3D) {
- sprintf(warn_buf,"With pbc=%s you should use ewald_geometry=%s",
+ sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
epbc_names[ir->ePBC],eewg_names[eewg3DC]);
warning(wi,warn_buf);
}
/* This check avoids extra pbc coding for exclusion corrections */
- sprintf(err_buf,"wall_ewald_zfac should be >= 2");
+ sprintf(err_buf,"wall-ewald-zfac should be >= 2");
CHECK(ir->wall_ewald_zfac < 2);
}
if (EVDW_SWITCHED(ir->vdwtype)) {
- sprintf(err_buf,"With vdwtype = %s rvdw_switch must be < rvdw",
+ sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw",
evdw_names[ir->vdwtype]);
CHECK(ir->rvdw_switch >= ir->rvdw);
} else if (ir->vdwtype == evdwCUT) {
ir->implicit_solvent=eisGBSA;
fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
"Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
- "setting implicit_solvent value to \"GBSA\" in input section.\n");
+ "setting implicit-solvent value to \"GBSA\" in input section.\n");
}
if(ir->sa_algorithm==esaSTILL)
}
}
+
+ if (ir->bAdress && !EI_SD(ir->eI)){
+ warning_error(wi,"AdresS simulation supports only stochastic dynamics");
+ }
+ if (ir->bAdress && ir->epc != epcNO){
+ warning_error(wi,"AdresS simulation does not support pressure coupling");
+ }
+ if (ir->bAdress && (EEL_PME(ir->coulombtype))){
+ warning_error(wi,"AdresS simulation does not support long-range electrostatics");
+ }
+
}
-static int str_nelem(const char *str,int maxptr,char *ptr[])
+int str_nelem(const char *str,int maxptr,char *ptr[])
{
int np=0;
char *copy0,*copy;
nstr = str_nelem(wall_density,MAXPTR,names);
if (nstr != ir->nwall)
{
- gmx_fatal(FARGS,"Expected %d elements for wall_density, found %d",ir->nwall,nstr);
+ gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
}
for(i=0; i<ir->nwall; i++)
{
sscanf(names[i],"%lf",&dbl);
if (dbl <= 0)
{
- gmx_fatal(FARGS,"wall_density[%d] = %f\n",i,dbl);
+ gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
}
ir->wall_density[i] = dbl;
}
RTYPE ("dt", ir->delta_t, 0.001);
STEPTYPE ("nsteps", ir->nsteps, 0);
CTYPE ("For exact run continuation or redoing part of a run");
- STEPTYPE ("init_step",ir->init_step, 0);
+ STEPTYPE ("init-step",ir->init_step, 0);
CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
- ITYPE ("simulation_part", ir->simulation_part, 1);
+ ITYPE ("simulation-part", ir->simulation_part, 1);
CTYPE ("mode for center of mass motion removal");
EETYPE("comm-mode", ir->comm_mode, ecm_names);
CTYPE ("number of steps for center of mass motion removal");
CTYPE ("Force tolerance and initial step-size");
RTYPE ("emtol", ir->em_tol, 10.0);
RTYPE ("emstep", ir->em_stepsize,0.01);
- CTYPE ("Max number of iterations in relax_shells");
+ CTYPE ("Max number of iterations in relax-shells");
ITYPE ("niter", ir->niter, 20);
CTYPE ("Step size (ps^2) for minimization of flexible constraints");
RTYPE ("fcstep", ir->fc_stepsize, 0);
ir->ndelta = 2;
CTYPE ("Periodic boundary conditions: xyz, no, xy");
EETYPE("pbc", ir->ePBC, epbc_names);
- EETYPE("periodic_molecules", ir->bPeriodicMols, yesno_names);
+ EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
CTYPE ("nblist cut-off");
RTYPE ("rlist", ir->rlist, 1.0);
CTYPE ("long-range cut-off for switched potentials");
RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
CTYPE ("Relative dielectric constant for the medium and the reaction field");
- RTYPE ("epsilon_r", ir->epsilon_r, 1.0);
- RTYPE ("epsilon_rf", ir->epsilon_rf, 1.0);
+ RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
+ RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
CTYPE ("Method for doing Van der Waals");
EETYPE("vdw-type", ir->vdwtype, evdw_names);
CTYPE ("cut-off lengths");
CTYPE ("Extension of the potential lookup tables beyond the cut-off");
RTYPE ("table-extension", ir->tabext, 1.0);
CTYPE ("Seperate tables between energy group pairs");
- STYPE ("energygrp_table", egptable, NULL);
+ STYPE ("energygrp-table", egptable, NULL);
CTYPE ("Spacing for the PME/PPPM FFT grid");
RTYPE ("fourierspacing", opts->fourierspacing,0.12);
CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
- ITYPE ("fourier_nx", ir->nkx, 0);
- ITYPE ("fourier_ny", ir->nky, 0);
- ITYPE ("fourier_nz", ir->nkz, 0);
+ ITYPE ("fourier-nx", ir->nkx, 0);
+ ITYPE ("fourier-ny", ir->nky, 0);
+ ITYPE ("fourier-nz", ir->nkz, 0);
CTYPE ("EWALD/PME/PPPM parameters");
- ITYPE ("pme_order", ir->pme_order, 4);
- RTYPE ("ewald_rtol", ir->ewald_rtol, 0.00001);
- EETYPE("ewald_geometry", ir->ewald_geometry, eewg_names);
- RTYPE ("epsilon_surface", ir->epsilon_surface, 0.0);
- EETYPE("optimize_fft",ir->bOptFFT, yesno_names);
+ ITYPE ("pme-order", ir->pme_order, 4);
+ RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
+ EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
+ RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
+ EETYPE("optimize-fft",ir->bOptFFT, yesno_names);
CCTYPE("IMPLICIT SOLVENT ALGORITHM");
- EETYPE("implicit_solvent", ir->implicit_solvent, eis_names);
+ EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
CTYPE ("Algorithm for calculating Born radii");
- EETYPE("gb_algorithm", ir->gb_algorithm, egb_names);
+ EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
CTYPE ("Frequency of calculating the Born radii inside rlist");
ITYPE ("nstgbradii", ir->nstgbradii, 1);
CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
CTYPE ("between rlist and rgbradii is updated every nstlist steps");
RTYPE ("rgbradii", ir->rgbradii, 1.0);
CTYPE ("Dielectric coefficient of the implicit solvent");
- RTYPE ("gb_epsilon_solvent",ir->gb_epsilon_solvent, 80.0);
+ RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
CTYPE ("Salt concentration in M for Generalized Born models");
- RTYPE ("gb_saltconc", ir->gb_saltconc, 0.0);
+ RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
- RTYPE ("gb_obc_alpha", ir->gb_obc_alpha, 1.0);
- RTYPE ("gb_obc_beta", ir->gb_obc_beta, 0.8);
- RTYPE ("gb_obc_gamma", ir->gb_obc_gamma, 4.85);
- RTYPE ("gb_dielectric_offset", ir->gb_dielectric_offset, 0.009);
- EETYPE("sa_algorithm", ir->sa_algorithm, esa_names);
+ RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
+ RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
+ RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
+ RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
+ EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
- RTYPE ("sa_surface_tension", ir->sa_surface_tension, -1);
+ RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
/* Coupling stuff */
CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
CTYPE ("Time constant (ps) and reference temperature (K)");
STYPE ("tau-t", tau_t, NULL);
STYPE ("ref-t", ref_t, NULL);
- CTYPE ("Pressure coupling");
- EETYPE("Pcoupl", ir->epc, epcoupl_names);
- EETYPE("Pcoupltype", ir->epct, epcoupltype_names);
+ CTYPE ("pressure coupling");
+ EETYPE("pcoupl", ir->epc, epcoupl_names);
+ EETYPE("pcoupltype", ir->epct, epcoupltype_names);
ITYPE ("nstpcouple", ir->nstpcouple, -1);
CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
RTYPE ("tau-p", ir->tau_p, 1.0);
STYPE ("compressibility", dumstr[0], NULL);
STYPE ("ref-p", dumstr[1], NULL);
CTYPE ("Scaling of reference coordinates, No, All or COM");
- EETYPE ("refcoord_scaling",ir->refcoord_scaling,erefscaling_names);
+ EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
CTYPE ("Random seed for Andersen thermostat");
- ITYPE ("andersen_seed", ir->andersen_seed, 815131);
+ ITYPE ("andersen-seed", ir->andersen_seed, 815131);
/* QMMM */
CCTYPE ("OPTIONS FOR QMMM calculations");
CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
STYPE ("annealing", anneal, NULL);
CTYPE ("Number of time points to use for specifying annealing in each group");
- STYPE ("annealing_npoints", anneal_npoints, NULL);
+ STYPE ("annealing-npoints", anneal_npoints, NULL);
CTYPE ("List of times at the annealing points for each group");
- STYPE ("annealing_time", anneal_time, NULL);
+ STYPE ("annealing-time", anneal_time, NULL);
CTYPE ("Temp. at each annealing point, for each group.");
- STYPE ("annealing_temp", anneal_temp, NULL);
+ STYPE ("annealing-temp", anneal_temp, NULL);
/* Startup run */
CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
/* Energy group exclusions */
CCTYPE ("ENERGY GROUP EXCLUSIONS");
CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
- STYPE ("energygrp_excl", egpexcl, NULL);
+ STYPE ("energygrp-excl", egpexcl, NULL);
/* Walls */
CCTYPE ("WALLS");
CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
ITYPE ("nwall", ir->nwall, 0);
- EETYPE("wall_type", ir->wall_type, ewt_names);
- RTYPE ("wall_r_linpot", ir->wall_r_linpot, -1);
- STYPE ("wall_atomtype", wall_atomtype, NULL);
- STYPE ("wall_density", wall_density, NULL);
- RTYPE ("wall_ewald_zfac", ir->wall_ewald_zfac, 3);
+ EETYPE("wall-type", ir->wall_type, ewt_names);
+ RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
+ STYPE ("wall-atomtype", wall_atomtype, NULL);
+ STYPE ("wall-density", wall_density, NULL);
+ RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
/* COM pulling */
CCTYPE("COM PULLING");
- CTYPE("Pull type: no, umbrella, constraint or constant_force");
+ CTYPE("Pull type: no, umbrella, constraint or constant-force");
EETYPE("pull", ir->ePull, epull_names);
if (ir->ePull != epullNO) {
snew(ir->pull,1);
pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
}
+
+ /* Enforced rotation */
+ CCTYPE("ENFORCED ROTATION");
+ CTYPE("Enforced rotation: No or Yes");
+ EETYPE("rotation", ir->bRot, yesno_names);
+ if (ir->bRot) {
+ snew(ir->rot,1);
+ rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
+ }
/* Refinement */
CCTYPE("NMR refinement stuff");
EETYPE("free-energy", ir->efep, efep_names);
RTYPE ("init-lambda", ir->init_lambda,0.0);
RTYPE ("delta-lambda",ir->delta_lambda,0.0);
- STYPE ("foreign_lambda", foreign_lambda, NULL);
+ STYPE ("foreign-lambda", foreign_lambda, NULL);
RTYPE ("sc-alpha",ir->sc_alpha,0.0);
ITYPE ("sc-power",ir->sc_power,0);
RTYPE ("sc-sigma",ir->sc_sigma,0.3);
EETYPE("separate-dhdl-file", ir->separate_dhdl_file,
separate_dhdl_file_names);
EETYPE("dhdl-derivatives", ir->dhdl_derivatives, dhdl_derivatives_names);
- ITYPE ("dh_hist_size", ir->dh_hist_size, 0);
- RTYPE ("dh_hist_spacing", ir->dh_hist_spacing, 0.1);
+ ITYPE ("dh-hist-size", ir->dh_hist_size, 0);
+ RTYPE ("dh-hist-spacing", ir->dh_hist_spacing, 0.1);
STYPE ("couple-moltype", couple_moltype, NULL);
EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
STYPE ("E-z", efield_z, NULL);
STYPE ("E-zt", efield_zt, NULL);
+ /* AdResS defined thingies */
+ CCTYPE ("AdResS parameters");
+ EETYPE("adress", ir->bAdress, yesno_names);
+ if (ir->bAdress) {
+ snew(ir->adress,1);
+ read_adressparams(&ninp,&inp,ir->adress,wi);
+ }
+
/* User defined thingies */
CCTYPE ("User defined thingies");
STYPE ("user1-grps", user1, NULL);
warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
}
} else {
- warning(wi,"Can not couple a molecule with free_energy = no");
+ warning(wi,"Can not couple a molecule with free-energy = no");
}
}
nref_t = str_nelem(ref_t,MAXPTR,ptr2);
ntcg = str_nelem(tcgrps,MAXPTR,ptr3);
if ((ntau_t != ntcg) || (nref_t != ntcg)) {
- gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref_t values and "
- "%d tau_t values",ntcg,nref_t,ntau_t);
+ gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
+ "%d tau-t values",ntcg,nref_t,ntau_t);
}
bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
snew(ir->opts.tau_t,nr);
snew(ir->opts.ref_t,nr);
if (ir->eI==eiBD && ir->bd_fric==0) {
- fprintf(stderr,"bd_fric=0, so tau_t will be used as the inverse friction constant(s)\n");
+ fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
}
if (bSetTCpar)
{
if (nr != nref_t)
{
- gmx_fatal(FARGS,"Not enough ref_t and tau_t values!");
+ gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
}
tau_min = 1e20;
ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
{
- sprintf(warn_buf,"With integrator %s tau_t should be larger than 0",ei_names[ir->eI]);
+ sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
warning_error(wi,warn_buf);
}
if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
{
if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
{
- sprintf(warn_buf,"For proper integration of the %s thermostat, tau_t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
+ sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
ETCOUPLTYPE(ir->etc),
tau_min,nstcmin,
ir->nsttcouple*ir->delta_t);
ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
if (ir->opts.ref_t[i] < 0)
{
- gmx_fatal(FARGS,"ref_t for group %d negative",i);
+ gmx_fatal(FARGS,"ref-t for group %d negative",i);
}
}
}
/* Read the other fields too */
nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
if(nSA_points!=nSA)
- gmx_fatal(FARGS,"Found %d annealing_npoints values for %d groups\n",nSA_points,nSA);
+ gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
for(k=0,i=0;i<nr;i++) {
ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
if(ir->opts.anneal_npoints[i]==1)
nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
if(nSA_time!=k)
- gmx_fatal(FARGS,"Found %d annealing_time values, wanter %d\n",nSA_time,k);
+ gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
if(nSA_temp!=k)
- gmx_fatal(FARGS,"Found %d annealing_temp values, wanted %d\n",nSA_temp,k);
+ gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
for(i=0,k=0;i<nr;i++) {
if (ir->ePull != epullNO) {
make_pull_groups(ir->pull,pull_grp,grps,gnames);
}
+
+ if (ir->bRot) {
+ make_rotation_groups(ir->rot,rot_grp,grps,gnames);
+ }
nacc = str_nelem(acc,MAXPTR,ptr1);
nacg = str_nelem(accgrps,MAXPTR,ptr2);
nr = groups->grps[egcENER].nr;
snew(ir->opts.egp_flags,nr*nr);
- bExcl = do_egp_flag(ir,groups,"energygrp_excl",egpexcl,EGP_EXCL);
+ bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
if (bExcl && EEL_FULL(ir->coulombtype))
warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
- bTable = do_egp_flag(ir,groups,"energygrp_table",egptable,EGP_TABLE);
+ bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
if (bTable && !(ir->vdwtype == evdwUSER) &&
!(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
!(ir->coulombtype == eelPMEUSERSWITCH))
decode_cos(efield_yt,&(ir->et[YY]),TRUE);
decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
-
+
+ if (ir->bAdress)
+ do_adress_index(ir->adress,groups,gnames,&(ir->opts),wi);
+
for(i=0; (i<grps->nr); i++)
sfree(gnames[i]);
sfree(gnames);
eel_names[eelGRF]);
CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
}
-
- if (ir->eI == eiSD1) {
- gdt_max = 0;
- for(i=0; (i<ir->opts.ngtc); i++)
- gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
- if (0.5*gdt_max > 0.0015) {
- sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta_t/tau_t = %g, you might want to switch to integrator %s\n",
- ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
- warning_note(wi,warn_buf);
- }
- }
+ if (ir->eI == eiSD1 &&
+ (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
+ gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
+ {
+ sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
+ warning_note(wi,warn_buf);
+ }
+
bAcc = FALSE;
for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
for(m=0; (m<DIM); m++) {
if (bConstr && ir->eConstrAlg == econtSHAKE) {
if (ir->shake_tol <= 0.0) {
- sprintf(warn_buf,"ERROR: shake_tol must be > 0 instead of %g\n",
+ sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
ir->shake_tol);
warning_error(wi,warn_buf);
}
}
if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
- sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs_order should be 8 or more.",ei_names[ir->eI]);
+ sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
warning_note(wi,warn_buf);
}
if (ir->epc==epcMTTK) {
#include "pbc.h"
#include "chargegroup.h"
#include "vec.h"
-#include "time.h"
+#include <time.h>
#include "nrnb.h"
#include "mshift.h"
#include "mdrun.h"
#include "trnio.h"
#include "xtcio.h"
#include "copyrite.h"
-
+#include "pull_rotation.h"
#include "mpelogging.h"
#include "domdec.h"
#include "partdec.h"
#ifdef GMX_LIB_MPI
#include <mpi.h>
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "tmpi.h"
#endif
+#include "adress.h"
#include "qmmm.h"
#if 0
double dt;
char buf[48];
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
if (!PAR(cr))
#endif
{
ir->delta_t/1000*24*60*60/runtime->time_per_step);
}
}
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
if (PAR(cr))
{
fprintf(out,"\n");
double mu[2*DIM];
gmx_bool bSepDVDL,bStateChanged,bNS,bFillGrid,bCalcCGCM,bBS;
gmx_bool bDoLongRange,bDoForces,bSepLRF;
+ gmx_bool bDoAdressWF;
matrix boxs;
real e,v,dvdl;
t_pbc pbc;
bDoLongRange = (fr->bTwinRange && bNS && (flags & GMX_FORCE_DOLR));
bDoForces = (flags & GMX_FORCE_FORCES);
bSepLRF = (bDoLongRange && bDoForces && (flags & GMX_FORCE_SEPLRF));
+ /* should probably move this to the forcerec since it doesn't change */
+ bDoAdressWF = ((fr->adress_type!=eAdressOff));
if (bStateChanged)
{
}
if (bStateChanged)
{
+
+ /* update adress weight beforehand */
+ if(bDoAdressWF)
+ {
+ /* need pbc for adress weight calculation with pbc_dx */
+ set_pbc(&pbc,inputrec->ePBC,box);
+ if(fr->adress_site == eAdressSITEcog)
+ {
+ update_adress_weights_cog(top->idef.iparams,top->idef.il,x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else if (fr->adress_site == eAdressSITEcom)
+ {
+ update_adress_weights_com(fplog,cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else if (fr->adress_site == eAdressSITEatomatom){
+ update_adress_weights_atom_per_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else
+ {
+ update_adress_weights_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ }
+
for(i=0; i<2; i++)
{
for(j=0;j<DIM;j++)
dd_force_flop_start(cr->dd,nrnb);
}
}
-
+
+ if (inputrec->bRot)
+ {
+ /* Enforced rotation has its own cycle counter that starts after the collective
+ * coordinates have been communicated. It is added to ddCyclF to allow
+ * for proper load-balancing */
+ wallcycle_start(wcycle,ewcROT);
+ do_rotation(cr,inputrec,box,x,t,step,wcycle,bNS);
+ wallcycle_stop(wcycle,ewcROT);
+ }
+
/* Start the force cycle counter.
* This counter is stopped in do_forcelow_level.
* No parallel communication should occur while this counter is running,
* since that will interfere with the dynamic load balancing.
*/
wallcycle_start(wcycle,ewcFORCE);
-
+
if (bDoForces)
{
/* Reset forces for which the virial is calculated separately:
if ((flags & GMX_FORCE_BONDED) && top->idef.il[F_POSRES].nr > 0)
{
- /* Position restraints always require full pbc */
- set_pbc(&pbc,inputrec->ePBC,box);
+ /* Position restraints always require full pbc. Check if we already did it for Adress */
+ if(!(bStateChanged && bDoAdressWF))
+ {
+ set_pbc(&pbc,inputrec->ePBC,box);
+ }
v = posres(top->idef.il[F_POSRES].nr,top->idef.il[F_POSRES].iatoms,
top->idef.iparams_posres,
(const rvec*)x,fr->f_novirsum,fr->vir_diag_posres,
start,homenr,mdatoms->chargeA,x,fr->f_novirsum,
inputrec->ex,inputrec->et,t);
}
+
+ if (bDoAdressWF && fr->adress_icor == eAdressICThermoForce)
+ {
+ /* Compute thermodynamic force in hybrid AdResS region */
+ adress_thermo_force(start,homenr,&(top->cgs),x,fr->f_novirsum,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
/* Communicate the forces */
if (PAR(cr))
if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
{
wallcycle_start(wcycle,ewcVSITESPREAD);
- spread_vsite_f(fplog,vsite,x,f,fr->fshift,nrnb,
+ spread_vsite_f(fplog,vsite,x,f,fr->fshift,FALSE,NULL,nrnb,
&top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
wallcycle_stop(wcycle,ewcVSITESPREAD);
if (bSepLRF)
{
wallcycle_start(wcycle,ewcVSITESPREAD);
- spread_vsite_f(fplog,vsite,x,fr->f_twin,NULL,
+ spread_vsite_f(fplog,vsite,x,fr->f_twin,NULL,FALSE,NULL,
nrnb,
&top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
wallcycle_stop(wcycle,ewcVSITESPREAD);
}
}
+ enerd->term[F_COM_PULL] = 0;
if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
{
/* Calculate the center of mass forces, this requires communication,
*/
set_pbc(&pbc,inputrec->ePBC,box);
dvdl = 0;
- enerd->term[F_COM_PULL] =
+ enerd->term[F_COM_PULL] +=
pull_potential(inputrec->ePull,inputrec->pull,mdatoms,&pbc,
cr,t,lambda,x,f,vir_force,&dvdl);
if (bSepDVDL)
}
enerd->dvdl_lin += dvdl;
}
+
+ /* Add the forces from enforced rotation potentials (if any) */
+ if (inputrec->bRot)
+ {
+ wallcycle_start(wcycle,ewcROTadd);
+ enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr,step,t);
+ wallcycle_stop(wcycle,ewcROTadd);
+ }
if (PAR(cr) && !(cr->duty & DUTY_PME))
{
* if the constructing atoms aren't local.
*/
wallcycle_start(wcycle,ewcVSITESPREAD);
- spread_vsite_f(fplog,vsite,x,fr->f_novirsum,NULL,nrnb,
+ spread_vsite_f(fplog,vsite,x,fr->f_novirsum,NULL,
+ (flags & GMX_FORCE_VIRIAL),fr->vir_el_recip,
+ nrnb,
&top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
wallcycle_stop(wcycle,ewcVSITESPREAD);
}
mtop,ir, (*outf)->fp_dhdl);
}
+ if (ir->bAdress)
+ {
+ please_cite(fplog,"Fritsch12");
+ please_cite(fplog,"Junghans10");
+ }
/* Initiate variables */
clear_mat(force_vir);
clear_mat(shake_vir);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff