check_function_exists(_fileno HAVE__FILENO)
check_function_exists(fileno HAVE_FILENO)
check_function_exists(_commit HAVE__COMMIT)
+check_function_exists(lstat HAVE_LSTAT)
check_function_exists(sigaction HAVE_SIGACTION)
include(CheckLibraryExists)
* The master node reads the file
* and communicates all the modified number of steps and the parallel setup,
* but not the state itself.
+ * When bAppend is set, lock the log file and truncate the existing output
+ * files so they can be appended.
+ * With bAppend and bForceAppend: truncate anyhow if the system does not
+ * support file locking.
*/
void load_checkpoint(const char *fn,FILE **fplog,
t_commrec *cr,gmx_bool bPartDecomp,ivec dd_nc,
t_inputrec *ir,t_state *state,gmx_bool *bReadRNG,
gmx_bool *bReadEkin,
- gmx_bool bTruncateOutputFiles);
+ gmx_bool bAppend,gmx_bool bForceAppend);
/* Read the state from checkpoint file.
* Arrays in state that are NULL are allocated.
#define MD_REPRODUCIBLE (1<<13)
#define MD_READ_RNG (1<<14)
#define MD_APPENDFILES (1<<15)
+#define MD_APPENDFILESSET (1<<21)
#define MD_KEEPANDNUMCPT (1<<16)
#define MD_READ_EKIN (1<<17)
#define MD_STARTFROMCPT (1<<18)
#endif
-#if (defined(WIN32) || defined( _WIN32 ) || defined(WIN64) || defined( _WIN64 ) || defined (__CYGWIN__))
+#if (defined(WIN32) || defined( _WIN32 ) || defined(WIN64) || defined( _WIN64 )) && !defined (__CYGWIN__) && !defined (__CYGWIN32__)
#define TMPI_NUMA_MALLOC
topology, this can be useful for normal mode analysis.</dd>
<dt>-DPOSRES</dt>
<dd>Will tell <tt>grompp</tt> to include posre.itp into your topology, used for
-<!--Idx-->position restraints<!--EIdx-->.</dd>
+<!--Idx-->position restraint<!--EIdx-->s.</dd>
</dl>
</dl>
<dt><b>integrator:</b> (Despite the name, this list includes algorithms that are not actually integrators. <tt>steep</tt> and all entries following it are in this category)</dt>
<dd><dl compact>
<dt><b>md</b></dt>
-<dd>A <!--Idx-->leap-frog<!--EIdx--> algorithm for integrating Newton's
-equations of motion.</dd>
+<dd>A leap-frog algorithm<!--QuietIdx-->leap-frog integrator<!--EQuietIdx-->
+for integrating Newton's equations of motion.</dd>
<dt><b>md-vv</b></dt>
<dd>A velocity Verlet algorithm for integrating Newton's equations of motion.
For constant NVE simulations started from corresponding points in the same trajectory, the trajectories
to the correction steps necessary it is not (yet) parallelized.
</dd>
<dt><b>nm</b></dt>
-<dd><!--Idx-->Normal mode analysis<!--EIdx--> is performed
+<dd>Normal mode analysis<!--QuietIdx-->normal-mode analysis<!--EQuietIdx--> is performed
on the structure in the <tt>tpr</tt> file. GROMACS should be
compiled in double precision.</dd>
<dt><b>tpi</b></dt>
<A NAME="em"><br>
<hr>
-<h3><!--Idx-->Energy minimization<!--EIdx--></h3>
+<h3>Energy minimization<!--QuietIdx-->energy minimization<!--EQuietIdx--></h3>
<dl>
<dt><b>emtol: (10.0) [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
<dd>the minimization is converged when the maximum force is smaller than
<A NAME="xmdrun"><br>
<hr>
-<h3><!--Idx-->Shell Molecular Dynamics<!--EIdx--></h3> When shells or
+<h3>Shell Molecular Dynamics<!--QuietIdx-->shell molecular dynamics<!--EQuietIdx--></h3>
+When shells or
flexible constraints are present in the system the positions of the shells
and the lengths of the flexible constraints are optimized at
every time step until either the RMS force on the shells and constraints
<A NAME="nl"><br>
<hr>
-<h3><!--Idx-->Neighbor searching<!--EIdx--></h3>
+<h3>Neighbor searching<!--QuietIdx-->neighbor searching<!--EQuietIdx--></h3>
<dl>
<dt><b>nstlist: (10) [steps]</b></dt>
<dd><dl compact>
<A NAME="el"><br>
<hr>
-<h3><!--Idx-->Electrostatics<!--EIdx--></h3>
+<h3>Electrostatics<!--QuietIdx-->electrostatics<!--EQuietIdx--></h3>
<dl>
<dt><b>coulombtype:</b></dt>
<dd><dl compact>
is identical to SPME, except that the influence function is optimized
for the grid. This gives a slight increase in accuracy.</dd>
-<dt><b><!--Idx-->Reaction-Field<!--EIdx--></b></dt>
+<dt><b>Reaction-Field electrostatics<!--QuietIdx-->reaction-field electrostatics<!--EQuietIdx--></b></dt>
<dd>Reaction field with Coulomb cut-off <b>rcoulomb</b>,
where <b>rcoulomb</b> ≥ <b>rlist</b>.
The dielectric constant beyond the cut-off is <b>epsilon-rf</b>.
<h3>Ewald</h3>
<dl>
<dt><b>fourierspacing: (0.12) [nm]</b></dt>
-<dd>The maximum grid spacing for the FFT grid when using PME or P3M.
-For ordinary Ewald the spacing times the box dimensions determines the
-highest magnitude to use in each direction. In all cases
-each direction can be overridden by entering a non-zero value for
-<b>fourier-n[xyz]</b>.
+<dd>For ordinary Ewald, the ratio of the box dimensions and the spacing
+determines a lower bound for the number of wave vectors to use in each
+(signed) direction. For PME and P3M, that ratio determines a lower bound
+for the number of Fourier-space grid points that will be used along that
+axis. In all cases, the number for each direction can be overridden by
+entering a non-zero value for <b>fourier_n[xyz]</b>.
For optimizing the relative load of the particle-particle interactions
-and the mesh part of PME it is useful to know that
+and the mesh part of PME, it is useful to know that
the accuracy of the electrostatics remains nearly constant
when the Coulomb cut-off and the PME grid spacing are scaled
by the same factor.</dd>
<A NAME="tc"><br>
<hr>
-<h3><!--Idx-->Temperature coupling<!--EIdx--></h3>
+<h3>Temperature coupling<!--QuietIdx-->temperature coupling<!--EQuietIdx--></h3>
<dl>
<dt><b>tcoupl:</b></dt>
<A NAME="pc"><br>
<hr>
-<h3><!--Idx-->Pressure coupling<!--EIdx--></h3>
+<h3>Pressure coupling<!--QuietIdx-->pressure coupling<!--EQuietIdx--></h3>
<dl>
<dt><b>pcoupl:</b></dt>
<A NAME="sa"><br>
<hr>
-<h3><!--Idx-->Simulated annealing<!--EIdx--></h3>
+<h3>Simulated annealing<!--QuietIdx-->simulated annealing<!--EQuietIdx--></h3>
Simulated annealing is controlled separately for each temperature group in GROMACS. The reference temperature is a piecewise linear function, but you can use an arbitrary number of points for each group, and choose either a single sequence or a periodic behaviour for each group. The actual annealing is performed by dynamically changing the reference temperature used in the thermostat algorithm selected, so remember that the system will usually not instantaneously reach the reference temperature!
<dl>
<h3>Bonds</h3>
<dl>
-<dt><b><!--Idx-->constraints<!--EIdx-->:</b></dt>
+<dt><b>constraints<!--QuietIdx-->constraint algorithms<!--QuietEIdx-->:</b></dt>
<dd><dl compact>
<dt><b>none</b></dt>
<dd>No constraints except for those defined explicitly in the topology,
<A NAME="walls"><br>
<hr>
-<h3><!--Idx-->Walls<!--EIdx--></h3>
+<h3>Walls<!--QuietIdx-->walls<!--EQuietIdx--></h3>
<dl>
<dt><b>nwall: 0</b></dt>
<dd>When set to <b>1</b> there is a wall at <tt>z=0</tt>, when set to <b>2</b>
<dt><b>disre:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
-<dd>no <!--Idx-->distance restraints<!--EIdx--> (ignore distance
-restraint information in topology file)</dd>
+<dd>ignore <!--Idx-->distance restraint<!--EIdx--> information in topology file</dd>
<dt><b>simple</b></dt>
-<dd>simple (per-molecule) distance restraints,
-ensemble averaging can be performed with <tt>mdrun -multi</tt>
-where the environment variable <tt>GMX_DISRE_ENSEMBLE_SIZE</tt> sets the number
-of systems within each ensemble (usually equal to the <tt>mdrun -multi</tt> value)</dd>
+<dd>simple (per-molecule) distance restraints.
<dt><b>ensemble</b></dt>
-<dd>distance restraints over an ensemble of molecules in one simulation box,
-should only be used for special cases, such as dimers
-(this option is not fuctional in the current version of GROMACS)</dd>
+<dd>distance restraints over an ensemble of molecules in one
+simulation box. Normally, one would perform ensemble averaging over
+multiple subsystems, each in a separate box, using <tt>mdrun -multi</tt>;s
+upply <tt>topol0.tpr</tt>, <tt>topol1.tpr</tt>, ... with different
+coordinates and/or velocities.
+The environment variable <tt>GMX_DISRE_ENSEMBLE_SIZE</tt> sets the number
+of systems within each ensemble (usually equal to the <tt>mdrun -multi</tt> value).</dd>
+</dd>
</dl></dd>
<dt><b>disre-weighting:</b></dt>
<dd><dl compact>
+<dt><b>equal</b> (default)</dt>
+<dd>divide the restraint force equally over all atom pairs in the restraint</dd>
<dt><b>conservative</b></dt>
<dd>the forces are the derivative of the restraint potential,
-this results in an r<sup>-7</sup> weighting of the atom pairs</dd>
-<dt><b>equal</b></dt>
-<dd>divide the restraint force equally over all atom pairs in the restraint</dd>
+this results in an r<sup>-7</sup> weighting of the atom pairs.
+The forces are conservative when <tt>disre-tau</tt> is zero.</dd>
</dl></dd>
<dt><b>disre-mixed:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>the violation used in the calculation of the restraint force is the
-time averaged violation </dd>
+time-averaged violation </dd>
<dt><b>yes</b></dt>
<dd>the violation used in the calculation of the restraint force is the
-square root of the time averaged violation times the instantaneous violation </dd>
+square root of the product of the time-averaged violation and the instantaneous violation</dd>
</dl></dd>
<dt><b>disre-fc: (1000) [kJ mol<sup>-1</sup> nm<sup>-2</sup>]</b></dt>
<dd>force constant for distance restraints, which is multiplied by a
-(possibly) different factor for each restraint</dd>
+(possibly) different factor for each restraint given in the <tt>fac</tt>
+column of the interaction in the topology file.</dd>
<dt><b>disre-tau: (0) [ps]</b></dt>
-<dd>time constant for distance restraints running average</dd>
+<dd>time constant for distance restraints running average. A value of zero turns off time averaging.</dd>
<dt><b>nstdisreout: (100) [steps]</b></dt>
-<dd>frequency to write the running time averaged and instantaneous distances
-of all atom pairs involved in restraints to the energy file
+<dd>period between steps when the running time-averaged and instantaneous distances
+of all atom pairs involved in restraints are written to the energy file
(can make the energy file very large)</dd>
<A NAME="nmr2">
<dt><b>orire:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
-<dd>no <!--Idx-->orientation restraints<!--EIdx--> (ignore orientation
-restraint information in topology file)</dd>
+<dd>ignore <!--Idx-->orientation restraint<!--EIdx--> information in topology file</dd>
<dt><b>yes</b></dt>
<dd>use orientation restraints, ensemble averaging can be performed
with <tt>mdrun -multi</tt></dd>
</dl>
<dt><b>orire-fc: (0) [kJ mol]</b></dt>
<dd>force constant for orientation restraints, which is multiplied by a
-(possibly) different factor for each restraint, can be set to zero to
+(possibly) different weight factor for each restraint, can be set to zero to
obtain the orientations from a free simulation</dd>
<dt><b>orire-tau: (0) [ps]</b></dt>
-<dd>time constant for orientation restraints running average</dd>
+<dd>time constant for orientation restraints running average. A value of zero turns off time averaging.</dd>
<dt><b>orire-fitgrp: </b></dt>
-<dd>fit group for orientation restraining, for a protein backbone is a good
+<dd>fit group for orientation restraining. This group of atoms is used
+to determine the rotation <b>R</b> of the system with respect to the
+reference orientation. The reference orientation is the starting
+conformation of the first subsystem. For a protein, backbone is a reasonable
choice</dd>
<dt><b>nstorireout: (100) [steps]</b></dt>
-<dd>frequency to write the running time averaged and instantaneous orientations
-for all restraints and the molecular order tensor to the energy file
+<dd>period between steps when the running time-averaged and instantaneous orientations
+for all restraints, and the molecular order tensor are written to the energy file
(can make the energy file very large)</dd>
</dl>
<A NAME="free"><br>
<hr>
-<h3><!--Idx-->Free energy calculations<!--EIdx--></h3>
+<h3>Free energy calculations<!--QuietIdx-->free energy calculations<!--EQuietIdx--></h3>
<dl>
<dt><b>free-energy:</b></dt>
<A NAME="neq"><br>
<hr>
-<h3><!--Idx-->Non-equilibrium MD<!--EIdx--></h3>
+<h3>Non-equilibrium MD<!--QuietIdx-->non-equilibrium MD<!--EQuietIdx--></h3>
<dl>
<dt><b>acc-grps: </b></dt>
<A NAME="ef"><br>
<hr>
-<h3><!--Idx-->Electric field<!--EIdx-->s</h3>
+<h3>Electric fields<!--QuietIdx-->electric field<!--EQuietIdx--></h3>
<dl>
<dt><b>E-x ; E-y ; E-z:</b></dt>
<hr>
<A NAME="qmmm"><br>
-<h3><!--Idx-->Mixed quantum/classical molecular dynamics<!--EIdx--></h3>
+<h3>Mixed quantum/classical molecular dynamics<!--QuietIdx>QM/MM<!--EQuietIdx--></h3>
<dl>
<dt><b>QMMM:</b></dt>
/* Define to 1 if you have the _fileno() function. */
#cmakedefine HAVE__FILENO
+/* Define to 1 if you have the lstat() function. */
+#cmakedefine HAVE_LSTAT
+
/* Define to 1 if you have the <string.h> header file. */
#cmakedefine HAVE_STRING_H
fprintf(debug,"Entries in %s: %d\n",ap->db,ap->nprop);
if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW)) {
- printf("\n");
- printf("WARNING: masses and atomic (Van der Waals) radii will be determined\n");
- printf(" based on residue and atom names. These numbers can deviate\n");
- printf(" from the correct mass and radius of the atom type.\n");
- printf("\n");
+ printf("\n"
+ "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
+ " based on residue and atom names, since they could not be\n"
+ " definitively assigned from the information in your input\n"
+ " files. These guessed numbers might deviate from the mass\n"
+ " and radius of the atom type. Please check the output\n"
+ " files if necessary.\n\n");
aps->bWarned = TRUE;
}
}
t_commrec *cr,gmx_bool bPartDecomp,ivec dd_nc,
int eIntegrator,gmx_large_int_t *step,double *t,
t_state *state,gmx_bool *bReadRNG,gmx_bool *bReadEkin,
- int *simulation_part,gmx_bool bAppendOutputFiles)
+ int *simulation_part,
+ gmx_bool bAppendOutputFiles,gmx_bool bForceAppend)
{
t_fileio *fp;
int i,j,rc;
* locking
*/
gmx_fatal(FARGS,"The first output file should always be the log "
- "file but instead is: %s", outputfiles[0].filename);
+ "file but instead is: %s. Cannot do appending because of this condition.", outputfiles[0].filename);
}
for(i=0;i<nfiles;i++)
{
/* lock log file */
if (i==0)
{
+ /* Note that there are systems where the lock operation
+ * will succeed, but a second process can also lock the file.
+ * We should probably try to detect this.
+ */
#ifndef GMX_NATIVE_WINDOWS
if (fcntl(fileno(gmx_fio_getfp(chksum_file)), F_SETLK, &fl)
==-1)
if (_locking(fileno(gmx_fio_getfp(chksum_file)), _LK_NBLCK, LONG_MAX)==-1)
#endif
{
- if (errno!=EACCES && errno!=EAGAIN)
+ if (errno == ENOSYS)
{
- gmx_fatal(FARGS,"Failed to lock: %s. %s.",
- outputfiles[i].filename, strerror(errno));
+ if (!bForceAppend)
+ {
+ gmx_fatal(FARGS,"File locking is not supported on this system. Use -noappend or specify -append explicitly to append anyhow.");
+ }
+ else
+ {
+ fprintf(stderr,"\nNOTE: File locking is not supported on this system, will not lock %s\n\n",outputfiles[i].filename);
+ if (fplog)
+ {
+ fprintf(fplog,"\nNOTE: File locking not supported on this system, will not lock %s\n\n",outputfiles[i].filename);
+ }
+ }
}
- else
+ else if (errno == EACCES || errno == EAGAIN)
{
gmx_fatal(FARGS,"Failed to lock: %s. Already running "
"simulation?", outputfiles[i].filename);
}
+ else
+ {
+ gmx_fatal(FARGS,"Failed to lock: %s. %s.",
+ outputfiles[i].filename, strerror(errno));
+ }
}
}
if (gmx_fio_get_file_md5(chksum_file,outputfiles[i].offset,
digest) != outputfiles[i].chksum_size) /*at the end of the call the file position is at the end of the file*/
{
- gmx_fatal(FARGS,"Can't read %d bytes of '%s' to compute checksum. The file has been replaced or its contents has been modified.",
+ gmx_fatal(FARGS,"Can't read %d bytes of '%s' to compute checksum. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
outputfiles[i].chksum_size,
outputfiles[i].filename);
}
}
fprintf(debug,"\n");
}
- gmx_fatal(FARGS,"Checksum wrong for '%s'. The file has been replaced or its contents has been modified.",
+ gmx_fatal(FARGS,"Checksum wrong for '%s'. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
outputfiles[i].filename);
}
#endif
#endif
if(rc!=0)
{
- gmx_fatal(FARGS,"Truncation of file %s failed.",outputfiles[i].filename);
+ gmx_fatal(FARGS,"Truncation of file %s failed. Cannot do appending because of this failure.",outputfiles[i].filename);
}
}
}
void load_checkpoint(const char *fn,FILE **fplog,
t_commrec *cr,gmx_bool bPartDecomp,ivec dd_nc,
t_inputrec *ir,t_state *state,
- gmx_bool *bReadRNG,gmx_bool *bReadEkin,gmx_bool bAppend)
+ gmx_bool *bReadRNG,gmx_bool *bReadEkin,
+ gmx_bool bAppend,gmx_bool bForceAppend)
{
gmx_large_int_t step;
double t;
read_checkpoint(fn,fplog,
cr,bPartDecomp,dd_nc,
ir->eI,&step,&t,state,bReadRNG,bReadEkin,
- &ir->simulation_part,bAppend);
+ &ir->simulation_part,bAppend,bForceAppend);
}
if (PAR(cr)) {
gmx_bcast(sizeof(cr->npmenodes),&cr->npmenodes,cr);
}
}
+static gmx_bool gmx_is_file(const char *fname)
+{
+ FILE *test;
+
+ if (fname == NULL)
+ return FALSE;
+ test=fopen(fname,"r");
+ if (test == NULL)
+ {
+ return FALSE;
+ }
+ else
+ {
+ fclose(test);
+ /*Windows doesn't allow fopen of directory - so we don't need to check this seperately */
+ #if (!((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__))
+ {
+ int status;
+ struct stat st_buf;
+ #ifdef HAVE_LSTAT
+ status = lstat (fname, &st_buf);
+ #else
+ status = stat (fname, &st_buf);
+ #endif
+ if (status != 0 || !S_ISREG(st_buf.st_mode))
+ {
+ return FALSE;
+ }
+ }
+ #endif
+ return TRUE;
+ }
+}
+
gmx_bool gmx_fexist_master(const char *fname, t_commrec *cr)
{
pdum=getcwd(system_path,sizeof(system_path)-1);
#endif
sprintf(full_path,"%s%c%s",system_path,DIR_SEPARATOR,bin_name);
- found = gmx_fexist(full_path);
+ found = gmx_is_file(full_path);
if (!found && (s=getenv("PATH")) != NULL)
{
char *dupped;
while(!found && (dir=gmx_strsep(&s, PATH_SEPARATOR)) != NULL)
{
sprintf(full_path,"%s%c%s",dir,DIR_SEPARATOR,bin_name);
- found = gmx_fexist(full_path);
+ found = gmx_is_file(full_path);
}
sfree(dupped);
}
np = interaction_function[ftype].nratoms;
if (ia[1] >= at_start && ia[1] < at_end) {
- if (ia[np] >= at_end)
+ if (ia[np] >= at_end || (ftype == F_SETTLE && ia[1]+2 >= at_end))
gmx_fatal(FARGS,
"Molecule in topology has atom numbers below and "
"above natoms (%d).\n"
{
fprintf(stderr,
"\n"
- "WARNING: if there are broken molecules in the trajectory file,\n"
- " they can not be made whole without a run input file\n\n");
+ "WARNING: If there are molecules in the input trajectory file\n"
+ " that are broken across periodic boundaries, they\n"
+ " cannot be made whole (or treated as whole) without\n"
+ " you providing a run input file.\n\n");
}
return gpbc;
#include "xvgr.h"
#include "gmxfio.h"
+#ifdef GMX_THREADS
#include "thread_mpi.h"
+#endif
#define p2(x) ((x)*(x))
#define p3(x) ((x)*(x)*(x))
}
else {
get_stx_coordnum(infile,&natoms);
- init_t_atoms(&top->atoms,natoms,FALSE);
+ init_t_atoms(&top->atoms,natoms,(fn2ftp(infile) == efPDB));
if (x == NULL)
{
snew(x,1);
/* I'm assuming we need global communication the first time! MRS */
cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
+ | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
| (bVV ? CGLO_PRESSURE:0)
| (bVV ? CGLO_CONSTRAINT:0)
| (bRerunMD ? CGLO_RERUNMD:0)
NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
constr,NULL,FALSE,state->box,
top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
- cglo_flags &~ CGLO_PRESSURE);
+ cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
}
/* Calculate the initial half step temperature, and save the ekinh_old */
/* these CGLO_ options remain the same throughout the iteration */
cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
- (bStopCM ? CGLO_STOPCM : 0) |
(bGStat ? CGLO_GSTAT : 0)
);
top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
cglo_flags
| CGLO_ENERGY
+ | (bStopCM ? CGLO_STOPCM : 0)
| (bTemp ? CGLO_TEMPERATURE:0)
| (bPres ? CGLO_PRESSURE : 0)
| (bPres ? CGLO_CONSTRAINT : 0)
top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
cglo_flags
| (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
| (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
| (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
| (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
Flags = Flags | (bRerunVSite ? MD_RERUN_VSITE : 0);
Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
Flags = Flags | (bAppendFiles ? MD_APPENDFILES : 0);
+ Flags = Flags | (opt2parg_bSet("-append", asize(pa),pa) ? MD_APPENDFILESSET : 0);
Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0);
Flags = Flags | (sim_part>1 ? MD_STARTFROMCPT : 0);
Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
cr,Flags & MD_PARTDEC,ddxyz,
inputrec,state,&bReadRNG,&bReadEkin,
- (Flags & MD_APPENDFILES));
+ (Flags & MD_APPENDFILES),
+ (Flags & MD_APPENDFILESSET));
if (bReadRNG)
{
real cellsize_limit,real cutoff_dd,
gmx_bool bInterCGBondeds,gmx_bool bInterCGMultiBody)
{
- int nnodes_div,ldiv;
+ gmx_large_int_t nnodes_div,ldiv;
real limit;
if (MASTER(cr))
}
debug_gmx();
-
- /* Calculate center of mass velocity if necessary, also parallellized */
- if (bStopCM && !bRerunMD && bEner)
- {
- calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
- state->x,state->v,vcm);
- }
}
- if (bTemp || bPres || bEner || bConstrain)
+ /* Calculate center of mass velocity if necessary, also parallellized */
+ if (bStopCM)
+ {
+ calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
+ state->x,state->v,vcm);
+ }
+
+ if (bTemp || bStopCM || bPres || bEner || bConstrain)
{
if (!bGStat)
{
wallcycle_start(wcycle,ewcMoveE);
GMX_MPE_LOG(ev_global_stat_start);
global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
- ir,ekind,constr,vcm,
+ ir,ekind,constr,bStopCM ? vcm : NULL,
gs != NULL ? eglsNR : 0,gs_buf,
top_global,state,
*bSumEkinhOld,flags);
mdatoms->massT,mdatoms->tmass,ekind->ekin);
}
- if (bEner) {
- /* Do center of mass motion removal */
- if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
- {
- check_cm_grp(fplog,vcm,ir,1);
- do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
- state->x,state->v,vcm);
- inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
- }
+ /* Do center of mass motion removal */
+ if (bStopCM)
+ {
+ check_cm_grp(fplog,vcm,ir,1);
+ do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
+ state->x,state->v,vcm);
+ inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
+ }
+ if (bEner)
+ {
/* Calculate the amplitude of the cosine velocity profile */
ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
}
break;
case F_POLARIZATION:
case F_ANHARM_POL:
- if (qS != atom[aS].qB)
- gmx_fatal(FARGS,"polarize can not be used with qA != qB");
+ if (!gmx_within_tol(qS, atom[aS].qB, GMX_REAL_EPS*10))
+ gmx_fatal(FARGS,"polarize can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS, atom[aS].qB, aS+1, mb+1);
shell[nsi].k += sqr(qS)*ONE_4PI_EPS0/
ffparams->iparams[type].polarize.alpha;
break;
case F_WATER_POL:
- if (qS != atom[aS].qB)
- gmx_fatal(FARGS,"water_pol can not be used with qA != qB");
+ if (!gmx_within_tol(qS, atom[aS].qB, GMX_REAL_EPS*10))
+ gmx_fatal(FARGS,"water_pol can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS, atom[aS].qB, aS+1, mb+1);
alpha = (ffparams->iparams[type].wpol.al_x+
ffparams->iparams[type].wpol.al_y+
ffparams->iparams[type].wpol.al_z)/3.0;
static void shell_pos_sd(FILE *log,rvec xcur[],rvec xnew[],rvec f[],
int ns,t_shell s[],int count)
{
+ const real step_scale_min = 0.8,
+ step_scale_increment = 0.2,
+ step_scale_max = 1.2,
+ step_scale_multiple = (step_scale_max - step_scale_min) / step_scale_increment;
int i,shell,d;
real dx,df,k_est;
#ifdef PRINT_STEP
for(d=0; d<DIM; d++) {
dx = xcur[shell][d] - s[i].xold[d];
df = f[shell][d] - s[i].fold[d];
- if (dx != 0 && df != 0) {
- k_est = -dx/df;
- if (k_est >= 2*s[i].step[d]) {
- s[i].step[d] *= 1.2;
- } else if (k_est <= 0) {
- s[i].step[d] *= 0.8;
- } else {
- s[i].step[d] = 0.8*s[i].step[d] + 0.2*k_est;
- }
- } else if (dx != 0) {
- s[i].step[d] *= 1.2;
- }
+ /* -dx/df gets used to generate an interpolated value, but would
+ * cause a NaN if df were binary-equal to zero. Values close to
+ * zero won't cause problems (because of the min() and max()), so
+ * just testing for binary inequality is OK. */
+ if (0.0 != df)
+ {
+ k_est = -dx/df;
+ /* Scale the step size by a factor interpolated from
+ * step_scale_min to step_scale_max, as k_est goes from 0 to
+ * step_scale_multiple * s[i].step[d] */
+ s[i].step[d] =
+ step_scale_min * s[i].step[d] +
+ step_scale_increment * min(step_scale_multiple * s[i].step[d], max(k_est, 0));
+ }
+ else
+ {
+ /* Here 0 == df */
+ if (gmx_numzero(dx)) /* 0 == dx */
+ {
+ /* Likely this will never happen, but if it does just
+ * don't scale the step. */
+ }
+ else /* 0 != dx */
+ {
+ s[i].step[d] *= step_scale_max;
+ }
+ }
#ifdef PRINT_STEP
step_min = min(step_min,s[i].step[d]);
step_max = max(step_max,s[i].step[d]);
{
int i,m,start,end;
gmx_large_int_t step;
- double mass,tmass,vcm[4];
real dt=ir->delta_t;
real dvdlambda;
rvec *savex;
}
}
- for(m=0; (m<4); m++)
- vcm[m] = 0;
- for(i=start; i<end; i++) {
- mass = md->massT[i];
- for(m=0; m<DIM; m++) {
- vcm[m] += state->v[i][m]*mass;
- }
- vcm[3] += mass;
- }
-
- if (ir->nstcomm != 0 || debug) {
- /* Compute the global sum of vcm */
- if (debug)
- fprintf(debug,"vcm: %8.3f %8.3f %8.3f,"
- " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],vcm[3]);
- if (PAR(cr))
- gmx_sumd(4,vcm,cr);
- tmass = vcm[3];
- for(m=0; (m<DIM); m++)
- vcm[m] /= tmass;
- if (debug)
- fprintf(debug,"vcm: %8.3f %8.3f %8.3f,"
- " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],tmass);
- if (ir->nstcomm != 0) {
- /* Now we have the velocity of center of mass, let's remove it */
- for(i=start; (i<end); i++) {
- for(m=0; (m<DIM); m++)
- state->v[i][m] -= vcm[m];
- }
-
- }
- }
sfree(savex);
}
}
}
}
+ }
- if (vcm)
+ if (vcm)
+ {
+ icm = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
+ where();
+ imass = add_binr(rb,vcm->nr,vcm->group_mass);
+ where();
+ if (vcm->mode == ecmANGULAR)
{
- icm = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
+ icj = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
where();
- imass = add_binr(rb,vcm->nr,vcm->group_mass);
+ icx = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
+ where();
+ ici = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
where();
- if (vcm->mode == ecmANGULAR)
- {
- icj = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
- where();
- icx = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
- where();
- ici = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
- where();
- }
}
}
+
if (DOMAINDECOMP(cr))
{
nb = cr->dd->nbonded_local;
extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
}
}
- /* should this be here, or with ekin?*/
- if (vcm)
- {
- extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
- where();
- extract_binr(rb,imass,vcm->nr,vcm->group_mass);
- where();
- if (vcm->mode == ecmANGULAR)
- {
- extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
- where();
- extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
- where();
- extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
- where();
- }
- }
if (DOMAINDECOMP(cr))
{
extract_bind(rb,inb,1,&nb);
}
}
+ if (vcm)
+ {
+ extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
+ where();
+ extract_binr(rb,imass,vcm->nr,vcm->group_mass);
+ where();
+ if (vcm->mode == ecmANGULAR)
+ {
+ extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
+ where();
+ extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
+ where();
+ extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
+ where();
+ }
+ }
+
if (nsig > 0)
{
extract_binr(rb,isig,nsig,sig);
"To convert a truncated octrahedron file produced by a package which uses",
"a cubic box with the corners cut off (such as GROMOS), use:[BR]",
"[TT]editconf -f in -rotate 0 45 35.264 -bt o -box veclen -o out[tt][BR]",
- "where [TT]veclen[tt] is the size of the cubic box times sqrt(3)/2." };
+ "where [TT]veclen[tt] is the size of the cubic box times [SQRT]3[sqrt]/2." };
const char *bugs[] =
{
"For complex molecules, the periodicity removal routine may break down, "
}
if (check_box(epbcXYZ,box))
- printf("\nWARNING: %s\n",check_box(epbcXYZ,box));
+ printf("\nWARNING: %s\n"
+ "See the GROMACS manual for a description of the requirements that\n"
+ "must be satisfied by descriptions of simulation cells.\n",
+ check_box(epbcXYZ,box));
if (bDist && btype[0][0]=='t')
{
load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
cr,Flags & MD_PARTDEC,ddxyz,
inputrec,state,&bReadRNG,&bReadEkin,
- (Flags & MD_APPENDFILES));
+ (Flags & MD_APPENDFILES),
+ (Flags & MD_APPENDFILESSET));
if (bReadRNG)
{
Flags = Flags | (bRerunVSite ? MD_RERUN_VSITE : 0);
Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
Flags = Flags | (bAppendFiles ? MD_APPENDFILES : 0);
+ Flags = Flags | (opt2parg_bSet("-append", asize(pa),pa) ? MD_APPENDFILESSET : 0);
Flags = Flags | (sim_part>1 ? MD_STARTFROMCPT : 0);
Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
"be combined with output from other programs and/or external",
"analysis programs to calculate more complex things.",
"Any combination of the output options is possible, but note",
- "that [TT]-om[tt] only operates on the first selection.[PAR]",
+ "that [TT]-om[tt] only operates on the first selection.",
+ "[TT]-os[tt] is the default output option if none is selected.[PAR]",
"With [TT]-os[tt], calculates the number of positions in each",
"selection for each frame. With [TT]-norm[tt], the output is",
"between 0 and 1 and describes the fraction from the maximum",
static real sbin=0,rmax=2,rbin=0.01,mmax=0,rint=0;
t_pargs pa[] = {
{ "-sqrt", FALSE, etREAL,{&sbin},
- "Use [SQRT]t[sqrt] on the matrix axis which binspacing # in sqrt(ps)" },
+ "Use [SQRT]t[sqrt] on the matrix axis which binspacing # in [SQRT]ps[sqrt]" },
{ "-fm", FALSE, etINT, {&fmmax},
"Number of frames in the matrix, 0 is plot all" },
{ "-rmax", FALSE, etREAL, {&rmax},
char ***atomname,t_atom atom[],
t_resinfo *resinfo)
{
- static const char * bb_nm[] = { "N", "H", "CA", "C", "O" };
+ static const char * bb_nm[] = { "N", "H", "CA", "C", "O", "HN" };
#define NBB asize(bb_nm)
t_bb *bb;
char *grpname;
bb[ri].N=ai;
break;
case 1:
+ case 5:
+ /* No attempt to address the case where some weird input has both H and HN atoms in the group */
bb[ri].H=ai;
break;
case 2:
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff