c42b1

stringm

String method for the study of conformational transitions V. Ovchinnikov (ovchinnv@georgetown.edu) * Description / Introduction to the string method suite * Invocation / String method invocation * SM0K / String method at zero temperature * SMCV / String method in collective variables * FTSM / Finite-temperature string method * Supporting files / Test cases and supporting files * References / References
Top Introduction. The string method module is an algorithm for finding paths between different configurations of a molecular system. It is a `chain-of-states' method, similar in principle to the Replica Path and Nudged Elastic Band (NEB) methods, in which a continuous transition path is `discretized' into a finite collection of system replicas under the constraint of approximately equal spacing using a predefined distance metric such as root-mean-square distance (RMSD)between adjacent replicas. Algorithm. The string methods are implemented according to the reference at the end of this document. For each replica, the 0-K string algorithm performs steepest descent (SD) evolution on the potential energy landscape (i.e. minimization) followed by a collective `reparametrization' to maintain equal spacing between adjacent replicas. The String method in collective variables (SMCV) by default performs SD minimization on a multidimensional _free energy_ landscape of the chosen collective variables. The required free energy gradients are samples by molecular dynamics with restraints. The finite-temperature string method (FTSM) samples hyperplanes locally perpendicular to an average path, computing free energies by integrating the average force acting on the hyperplanes, or by sampling a tessellation of the conformational space defined on the basis of the average path. Implementation notes The string method is a parallel simulation technique and requires that message-passing interface and the keyword 'STRINGM' to compile the string method source code): $> ./install.com gnu gfortran stringm Each bead on the string is represented by a group of processors which communicate internally via the standard communicator MPI_COMM_LOCAL (which is the same as COMM_CHARMM inside the program). The exchange of information between each processor group, as required during, e.g., the reparametrization step, occurs via an additional string communicator (MPI_COMM_STRNG), whose nodelist consists of the root processor in each local group. This communication scheme is set up using the module multicom at the beginning of each string calculation (see command references below). Limitations The zero-temperature string method (SM0K) is designed for the steepset descent minimizer, although other minimizers can be used 'manually' by explicitly invoking reparametrization and minimization alternately. The string methods that use dynamics (SMCV and FTSM) are designed for the standard leapfrog integrator under the original by-atom parallelization scheme. Thus, other integrators and parallel schemes (e.g. velocity verlet, DOMDEC, constant pressure integration) are not supported, even if the CHARMM code may compile with many different options. This retains the convenience of using a single executable for different types of simulations, but the compatibility restrictions should be kept in mind. Extensions of select parts of the string methods to other integrators/parallelizations may be implemented in the future.
Top Invocation of the string method modules. Introduction. The string method suite relies on a special communication scheme that is provided by the module multicom (MCOM). This module is to be called prior to any string calculation to define the message-passing communicators for communication within and between string nodes. (This implies that the string method requires the message passing interface, which in a CHARMM compilation is provided by the keyword M). A special command to MCOM will create a series of local communicators to handle communication between different cpus associated with the same string image, and a communicator for communication between different string images. Subsequent calls to string method initialization routines will check that the appropriate communicators have been set up and print a status message or a warning. Syntax. MCOM STRIng <int> BY <int> STRIng [ ZERO <character*> | SM0K <character*> | COLV <character*> | FTSM <character*> | OPEN [UNIT <int>] - [NAME <character*>] - [{UNFOrmatted|FILE|FORMatted|CARD}] - [{READ|WRITe|APPEnd}] - CLOSe] Command description. MCOM STRIng <int>(i) BY <int>(j) Create MPI communicators for (j) string replicas, with (i) CPU cores assigned to each replica. The total number of CPUs in the parallel run (e.g. the number given to the mpirun command must be no less than (i)*(j), otherwise an error message will be printed. ------------------------------------------------------------------------ [SM0K|ZERO] Invoke the command parser for the zero-temperature string method (SM0K subcommands are documented elsewhere in this file) [SMCV|COLV] Invoke the command parser for the string method in collective variables (COLV subcommands are documented elsewhere) [FTSM] Invoke the zero-temperature string command parser (FTSM subcommands are documented elsewhere) [OPEN] Open a file on the root node of each string image. I/O specification is the same as that described for general CHARMM I/O (see io.doc). This STRIng OPEN command is different from the OPEN command because the latter will only open a file on the root note of the parallel run. Known issues : Execution of MCOM STRIng may require modification to sequence definitions (i.e. READ SEQUence). More precisely, e.g., the following snippet in the main script will not read a new sequence on all string replicas: read sequence card * new sequence [n] ALA ... The sequence will be read only on the global root node. This is because the string roots will only read a file if it has been opened with "STRIng OPEN", which is impossible to do above. To solve this issue, two methods are recommended: (1) read sequence prior to calling MCOM (2) put the sequence in a text file and open the text file with "STRIng OPEN", i.e. string open unit 50 seq.dat read sequ card unit 50 [CLOSe] Close file unit on the root node of each string image.
Top Zero-temperature (0K) String Method STRIng ZERO [INITialize]} | [REPArametrize [{DEFI <real>}] - [{ITER <int>}] - [{LINE | SPLI | BSPL | DST [{WNCT <real>}] | LIN2}] - [{ORIE <atom-selection>} [{MASS}] ] - [{MOVE <atom-selection>}] ] | [STATistics [{COUNt <int>]} - [{ENER <energy_terms> [ENAM <character*>]}] - [{RMSD [RNAM <character*>] [{RAPP}]} ] - [{RMSA [RANM <character*>] [{RAAP}]} ] - [{DELS [DNAM <character*>] [{DAPP}]} ] - [{ORIE}] - [{MASS}] - [{ARCL [ANAM <character*>] [{AAPP}]} ] - [{CURV [CVNM <character*>] [{CAPP}]} ] ] | [CONFormational Consistency}] | [CHIRality}] | [MINImize [{REPF <int>}] - [{STAF <int>}] - [{CONFF <int>}] - [{CHIRF <int>}] - [{NOPR}]] | [INTErpolate [{ METH [LINE|SPLI|BSPL]}] - [NIN <int>] - [NOUT <int>] - [{CRIN <character*>} | {DCDIN <character*> [{BEGIn <int>}] [{STEP <int>}] } ] - [OUTLIST <character*> | OUTBASE <character*> [{OUTI <int>}]] - [{PDB [{RESI}] | FILE | UNFO | CARD | FORM}] - [{ORIE <atom selection> [{MASS}]}] ] -------------------------------------------------------------------------- [INIT] Initialize the zero-temperature string module. Substitution parameters ?MESTRING and ?NSTRING will be set to correspond to the string rank of each CPU group, and to the total number of string replicas, respectively. -------------------------------------------------------------------------- [REPA] Call reparametrization (Repa) module. [LINE] Perform simple linear interpolation. This is the most stable and preferred method. [SPLI] Perform cubic spline interpolation. This method has very low dissipation, and may cause the string curves corresponding to very complex systems (and especially poorly resolved) systems to lengthen indefinitely, preventing convergence. [BSPL] Perform B-spline interpolation. This method introduces more dissipation than the LINE method, and produces a smoother string curve. [DST {WNCT <real> w \in [0..1] } ] Perform interpolation using the Discrete Sine Transform described in reference [3]. This method presently appears to offer no advantages over the other interpolation methods above, is not often used, and thus should be considered experimental. [WNCT <real>] specifies the wave number cut-off for the DST method. A value w implies that only the lowest w-fraction of wavenumbers will be used to reconstruct the path by the inverse sine transform. [LIN2] Exact (i.e. noniterative) linear interpolation in which Newton-Raphson iterations are used to redistribute images such length increments between adjacent images are exactly equal. Because the NR iterations may fail to converge for very jagged strings, the method is not generally recommended for long simulations. [DEFI <real>] maximum allowed RMSD error in the distance between adjacent replicas, defined as max_i RMSD[i,i+1]/RMSD[i,i-1]. The default value is 1.1. A value <= 1.0 will cause the reparametrization to proceed to the maximum number of iterations (see below). [ITER <int> ] maximum number of reparametrization iterations. The default value is 10. NOTE: iterative reparametrization will continue until the DEFI requirement is satisfied, or the maximum number of iterations is exceeded. 'ITER 0' is allowed, and will cause the code to compute information about the string parametrization without modifying it. [ORIE {SELE <atom-selection> END}] atom selection that indicates that the adjacent string replicas are to be superposed to miminize the RMSD between atoms in the atom selection. This option should typicaly be present for rigid-body-invariant systems. The [MASS] specifies that the orientation is to use mass-weighting. Prior to reparametrization, replica (i>1) will be rotated/translated such that the RMSD between the orientation sets of atoms of replicas (i) and (i-1) is minimized. Replica coordinates will be restored to the original frame of reference after reparametrization. [MOVE {SELE <atom-selection> END}] atom selection that specifies the parametrization of the string, i.e. those atoms to which reparametrization will be applied. Coordinates not included in teh selection will not be affected by reparametrization. The default option is to MOVE all atom coordinates. -------------------------------------------------------------------------- [STAT] Call statistics (Stat) module. Specify statistics options when arguments are present. When called without arguments, an instance of statistics output will be written out. [COUNt <int>] specifies iteration counted for statistics output. The default value is 0. This number is incremented after every call to statistics, and is printed in the output files corresponding to RMSD, DELS, ARCL, and CURV. It is also appended to the base file name specified by the [ENAM <character*>] [ENER {<energy_terms>} [ENAM <character*>]] sets up output of energy values along the path. At each statistics call, the energy values corresponding to the terms in {<energy_terms>} will be listed in one file per iteration, with one line of output per replica. The list {<energy_terms>} can contain one or more energy substitution keywords described in energy.doc (Energy ). The current iteration will be appended to the base energy file name specified with [ENAM <character*>], followed by extension '.DAT'. If ENAM is omitted, energy output is directed to the output stream. This subcommand must be terminated by END. [RMSD RNAM <character*> [RAPP]] sets up output of RMSD values between the replica coordinates at the current iteration and the coordinates present in the comparison set (usually the initial string). At each call to string statistics, a line will be added to the file with the name specified in [RNAM <character*>], with the columns in the file corresponding to different replicas. If RAPP is specified, output will be appended to the file. If RNAM is omitted, RMSD output is directed to the output stream. This subcommand is useful for gauging convergence of the string. [RMSA RANM <character*> [RAAP]] sets up output of RMSD values between the replica coordinates at the current iteration and the time-averaged simulation coordinates of the replica. At each call to string statistics, a line will be added to the file with the name specified in [RNAM <character*>]. If RAAP is specified, output will be appended to the file. [DELS DNAM <character*> [DAPP]] sets up output of RMSD values between the replica coordinates at the current iteration and those at the previous iteration. At each call to string statistics, a line will be added to the file with the name specified with [DNAM <character*>]. If DAPP is specified, output will be appended to the file. [ORIE {MASS}] If ORIE is present, structures will be superposed to minimize the RMSD between simulation and reference, averaged, and previous structures. If MASS is present, orientation will use mass-weighting. These keywords affect the behavior of options RMSD, RMSA, and DELS only. [ARCL ANAM <character*> [AAPP]] sets up output of the distance between the adjacent replicas (such that their sum yields the string length) at the current iteration. At each call to string statistics, a line will be added to the file with the name specified with [ANAM <character*>]. If AAPP is specified, output will be appended to the file. [CURV CVNM <character*> [CAPP]] sets up output of the approximate string curvature at each string image. At each call to string statistics, a line will be added to the file with the name specified with [CVNM <character*>]. If CAPP is specified, output will be appended to the file. -------------------------------------------------------------------------- [MINI] Perform string energy minimization using steepest descent (other minimizers are supported indirectly at script level). [REPF <int>] Number of steps to skip between consecutive reparametrizations. REPA options must be set prior to calling MINI. [CHIRF <int>] Number of steps between consecutive calls to conformational consistency. CHIR options must be set prior to calling MINI. [CONFF <int> ] Number of steps between consecutive calls to CONFCons module. CONF options must be set prior to calling MINI. [STAF <int>] Number of reparametrization steps between consecutive calls to statistics module, i.e. the number of minimization steps between consecutive statistics to is <iSTATF> * <iREPF> . [NOPR] Flag to limit non-essential string output from various modules during minimization. -------------------------------------------------------------------------- [INTE] Interpolate a string onto a new string with different resolution. Useful for creating initial paths between known endpoints. [METH [LINE|SPLI|BSPL] ] interpolation method. Presently valid choices are LINE, SPLI, BSPL (described above for REPA). The default is LINE. [NIN] Number of coordinate sets in the existing string. [NOUT] Number of coordinate sets in the interpolated string. [CRIN <character*>] Name of text file which contains the filenames in which existing coordinates are stored, with one file name per line. This option is mutually exclusive with DCDIN. [DCDIN <character*> [{BEGIn <int>}] [{STEP <int>}] } ] Name of trajectory file which contains the existing string coordinates. BEGIN <int> specifies the frame number corresponding to the first string image (default: 1). STEP <int> specifies the stride between adjacent string images in the trajectory (default: 1). [OUTLIST <character*> ] Indicates that output file names corresponding the interpolated string are to be found in a text file with the provided name (one file name per line). This option is mutually exclusive with OUTBASE. [OUTBASE] Indicates that output file names corresponding to the interpolated string are to be of the form <outbase><i>.cor; i \in [<outi> .. <outi> + <nout> - 1], where : [OUTI <int>outi ] specifies the initial index (default: 0) The following options are passed to the CHARMM coordinate input/output routines: [PDB [{RESI}] ] Indicates that input/output coordinate files are in PDB format. RESI will pass the 'RESID' option to the coordinate reader. [FILE | UNFO ] Indicates that input/output coordinate files are UNFORMATTED. [CARD | FORM}] Indicates that input/output coordinate files are FORMATTED. [ORIE {MASS}] If ORIE is present, structures will be superposed to minimize the RMSD between adjacent string images prior to interpolation. If MASS is present, orientation will use mass-weighting. -------------------------------------------------------------------------- [CONF] Call conformational consistency (ConfCons) module. CONFcons options can be set prior to calling using the CONFCons interface described in corman.doc. -------------------------------------------------------------------------- [CHIR] Call chirality (Chir) module. CHIRality options can be set prior to calling using the CHIRality interface described elsewhere in this file. -------------------------------------------------------------------------- Chirality module Description. A simple chirality checker to find and correct chirality errors acording to prescribed rules, applied to an optional atom selection. The functionality is provided as part of the string module along with the conformational consistency (CONFcons) module to correct coordinate errors that occur when generating minimum energy paths using the 0K string method. Syntax. COOR CHIRality [{INIT}] [{FINA}] [{CHECK}] [{FIX}] [{NOFX}] [{NOPR}] [{PRIN}] [{HELP}] [{RULE}] [{PSEU}] [{PROP}] [{READ {UNIT} {ADD}}] {<atom selection>} Command reference. [INIT ] initialize Chirality (also, reinitialize, erasing old rules). [FINA ] finalize Chirality, deallocating all data. [CHECK] check residues specified in the atom selection for errors [FIX ] attempt to correct Chirality errors found during checking. [NOFX ] check for errors, but do not attempt to correct them [NOPR ] produce minimal output (useful during 0K string minimizations). [PRIN ] produce normal (verbose) output [HELP ] print succinct command reference [RULE ] print currently defined rules that will be applied to each matching residue. The rule syntax is decsribed below. [PSEU ] also check for pseudo-chirality errors (i.e. ordering of hydrogens in a methyl group). [PROP ] also check for pseudo-chirality in the Proline patch. [READ {UNIT<int>} {ADD} ] read rules from unit (or from the input stream if unit is omitted). ADD will append new rules to the existing set of rules. Atom selection is optional, and all atoms will be selected by default. When present, each residue with at least one atom contained in the atom selection will be checked (i.e. BYRES keyword is implied in the selection) Description of chirality rules. Each chirality rule is a character string that contains the following in order: Residue name <char*>, four atom types <char*>, a dihedral angle threshold for determining whether the chirality is incorrect <real>, an flag that determines the correct dihedral range, and names of two atoms which specify the bond to be inverted in case of a chirality error. For example, the rule 'ARG N C CB HA 0. -1 HA CA' indicates that the dihedral angle N-C-CB-HA in residue ARG must be less than zero Precisely, the inequality that is evaluated is [ phi(x) * (-1) > phi0 * (-1) ], where phi0=0 and the flag "-1" implicitly specifies the valid range, i.e. valid `side` of the inequality. If this rule is violated, then the position of atom HA will be reflected through the plane perpendicular to the bond CA-HA at CA. Known issues. For some pathological geometries, in particular those that are far from tetrahedral, the chirality module will not be able to identify the proper reflection. More precisely, this will occur if reflecting the corresponding atom in the rule (e.g. HA in the above rule) will not lead to a sign change in the dihedral. This could happen, e.g., if the reflection atom is located too close to its bonded neighbor (e.g. CA in the above rule). It is therefore strongly recommended to use the chirality module with geometries that are close to equilibrium. -------------------------------------------------------------------------- Conformational consistency module. Description. The CONFormational CONSistency module is applicable to groups of atoms that have rotational symmetry because some of the constituent atoms are indistinguishable, e.g. H1, H2, and H3 in a methyl group. In specifying the coordinates of a molecule with such a group, the atoms can be arranged in three completely equivalent ways, relative to a reference conformation. These are H1/H2/H3, H2/H3/H1, H3/H1/H2. They are related by rotations around the bond that involves the parent atom (e.g. methyl carbon). Note that the conformation H1/H3/H2 is not allowed because it corresponds to a change of (pseudo) chirality, which is assumed to be fixed (see chirality module) The purpose of this module is to identify all such symmetry groups in a molecular coordinate set and choose the rotation that brings the coordinates into the closest proximity (in the sense of a certain RMSD ) to the coordinates in a reference set. This is done according to a set of rules, which can be specified by the user; a default set for the charmm22 protein topology is defined in this module. A description of the algorithm and the rules is provided in corman.doc. When used with the 0K string method, the CONFcons module check consistency for each pair of adjacent structures along the string. -------------------------------------------------------------------------- Test cases. c40test/sm0k.inp c40test/confcons.inp -------------------------------------------------------------------------- SM0K References [1] W. E, W. Ren, and E. Vanden-Eijnden, J. Chem. Phys. 126, 164103 (2007). [2] V. Ovchinnikov, M. Karplus, and E. Vanden-Eijnden, J. Chem. Phys. 134, 085103 (2011). [3] I. Khavrutskii, K. Arora, and C. Brooks III, J. Phys. Chem. 125, 174108 (2006). If you find the SM0K code useful in preparing a manuscript, please cite the implementation reference 2 above.
Top String Method in Collective Variables (SMCV). ------------------------------------------------------------------------------------------------------------- Command Syntax. STRIng COLVar [INITialize [{MAXCV}] | [DONE] | [ADD [POSI_COM_X[{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}][FRAMe<int>] <atom_selection>]| [POSI_COM_Y[{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}][FRAMe<int>] <atom_selection>]| [POSI_COM_Z[{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}][FRAMe<int>] <atom_selection>]| [DIHE_COM [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}]<atom_selection><x4> ] | [ANGLE_COM [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}]<atom_selection><x3> ] | [DIST_COM [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}]<atom_selection><x2> ] | [ANGLVEC [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}] - [P1 [<atom_selection> | 3<real>]] [P2 [<atom_selection> | 3<real>]] - [P3 [<atom_selection> | 3<real>]] [P4 [<atom_selection> | 3<real>]] - [{FR1<int>}] [{FR2<int>}] ] | [FRAME <atom_selection>] | [QUATERNION [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}][{FRA1<int>}][{FRA2<int>}]] | [RMSD [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}] <atom_selection> - [{ORIE <atom_selection>}] [{MASS}] [{MAIN|COMP}]] | [DRMSD [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}] <atom_selection> - [{ORIE <atom_selection>}] [{MASS}] [{MAIN|COMP}]] | [PROJ [{FORCe<real>}][{<WEIGht<real>}][{GAMMa<real>}] <atom_selection> - [{ORIE <atom_selection>}] [{MASS}] [{MAIN|COMP}]]] | [LIST] | [CLEAr] | [REPArametrize [{DEFI <real>}] - [{ITER <int>}] - [{LINE | SPLI | BSPL | DST [{WNCT <real>}] | LIN2}] | [TEST [GRAD {STEP <real>} | PARA | MINV] ] | [FILL {MAIN|REF|OLD} {COMP} ] | [SWAP 2x[MAIN|REF|OLD] ] | [COPY 2x[MAIN|REF|OLD] ] | [PRINt[{ALL}] {UNIT <int>} {NAME <character*>} {COL}] | [READ [{ALL | SCOL<int> TCOL<int>}] {UNIT <int>} {NAME <character*>} {COL}] | [SET [{IND<int> | ALL }] - [{FORCe<real>}] [{GAMMa<real>}] [{WEIGht<real>}] [{ZVAL<real> [REP<int>] [{COL<int>}]}] ] | [VOROnoi [VCUT <real> [{REP<int>}]] | [VMAP [CALC] | - [PRINt[{UNIT <int>}] [NAME <character*>] ] - [READ [{UNIT <int>}] [NAME <character*>] ] - [CLEA] ] | - [READ [UNIT <int>] [NAME <character*>] ] - [PRINt[UNIT <int>] [NAME <character*>] ] ] | [FRAMe[CLEAr] | [RESEt] | [FILL [{COMP}] [{ALIGn}]] | [PRINt[{ALL}] [{UNIT<int>}] [{NAME<character*>}]] | [READ [{ALL}] [{UNIT<int>}] [{NAME<character*>}]] | [LIST] | [ALIGn [{RMSD | VORO}]] ] | [PARA [QUAT [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] |- [FRAM [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] |- [COLV [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] |- [MMAT [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] |- [VORO [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] ] | [MINV [LU|lu|DIAG|diag]] | [MMAT [CALC | - [PRINt[{INV}] [{UNIT<int>}] [NAME<character*>]] | - [READ [{INV}] [{UNIT<int>}] [NAME<character*>]] | - [FAST [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]] ] | [WEIGht[CALC| [PRINt[{UNIT<int>}] [NAME<character*>]] | - [READ [{UNIT<int>}] [NAME<character*>]] ] | [INTE [{INTERPCV [METH [LINE|SPLI|BSPL|LIN2]]}] - [NIN<int>] [CVIN<character*>] - [NOUT<int>][CVOUT<character*>] - [{COOR [CRIN<character*>][CROUT<character*>] - [{PDB [{RESI}] | FILE | UNFO | CARD | FORM}]]] | [DYNAmics[{STAF <int>}] - [{NOPR}] - [{VORO}] - [{REPA [REPF <int>]} ] - [{STAT [STAF]}] - [{HISF <int>] - [{EVOL [EVOF <int>] - [{EVOS <int>}] - [{EVOStep <real>}] - [{EXPO {MEMO <real>} | AVER {NAVE <int> }]] [{RSTR [REEQ <int>] {SMD|STEER}}] - [{REX [REXF <int>] [REXT <real>]}] - <regular dynamics options>] | [STATistics[{COUNt <int>]} - [{RMSD [RNAM <character*>] [{RAPP}]} ] - [{RMSA [RANM <character*>] [{RAAP}]} ] - [{DELS [DNAM <character*>] [{DAPP}]} ] - [{ARCL [ANAM <character*>] [{AAPP}]} ] - [{CURV [CVNM <character*>] [{CAPP}]} ] - [{FREE [FENM <character*>] [{FAPP}]} ] - [{COLV [CNAM <character*>] [{CVAP}] }] [{VORO [VNAM <character*>]} ] - [{VMAP [{VNAM <character*>}]} ] - [{VLOG [{VNAM <character*>}]} [{VOFF <int>}] [{VLAP}]}] - [{REXM [RXNM <character*>] [{RXOL <character*>}]}] - [{REXL [{RXNM <character*>}] [{ROFF <int>}] [{RXAP}]}] - [{FORC [FCNM <character*>] [{FCAP}] }] - [{MMAT [MNAM <character*>]}]] Command Description. -------------------------------------------------------------------------- [INIT] Initialize the string method in collective variables (SMCV). Substitution parameters ?MESTRING and ?NSTRING will be set to correspond to the string rank of each CPU group, and to the total number of string replicas, respectively. -------------------------------------------------------------------------- [DONE] Finalize SMCV module -------------------------------------------------------------------------- [ADD [<CV_specification>]] Define a collective variable. The following CV are supported. 1)POSI_COM_X : X-component of COM of atom group 2)POSI_COM_Y : Y-component of COM of atom group 3)POSI_COM_Z : Z-component of COM of atom group 4)DIST : Distance between two COM groups 5)ANGL : Angle between two COM groups 6)ANGLVEC : Angle between two vectors 7)DIHE : Dihedral angle between two COM groups 8)RMSD : Root-mean-square distance to a reference coordinate set of an atom group 9)DRMSD : Difference between RMSDs to two coordinate sets 10)PROJ : Fraction of distance between two coordiante sets (projected onto vector connecting the two sets, accounting for rotation) 11)QUATERNION : quaternion representation of rotation transformation between two frames. When a quaternion CV is specified, four CVs will actually be created, corresponding to the components of the quaternion. A new coordinate reference frame can be specified on the basis of an atom selection using ADD FRAME <atom_selection>; the frame axes will be defined as the eigenvectors of the moment of inertia tensor of the selected coordinates. The following (optional) keywords are valid for all CV: [{FORCe<real>}] Set harmonic force constant. [{<WEIGht<real>}] Set CV weight in distance computation (this is akin to atomic mass in the case of absolute position CV. [{GAMMa<real>}] Set the CV friction coefficient for the evolution of the CV in inverse AKMA time units (see EVOL keyword for the evolution equation) NOTE: Each CV definition requires a corresponding number of atom selections as indicated in teh syntax section above. RMSD, DRMSD and PROJ CVs take an optional orientation set specified by the keyword ORIE after the main (RMSD) selection. The MASS keyword causes mass-weighting to be used in the RMSD computation ans in best-fit superpositions. POSI_COM_[X|Y|Z] CVs take a optional frame number via `FRAMe <int>` (default 0 for absolute frame) in which the positions will be expressed. To use a non-absolute frame, such a frame must first be defined via `ADD FRAMe <atom_selection>`; frame indices correspond to the order in which frames have been defined. The ANGLVEC CV specifies the angle between two vectors. The vectors are defined using four coordinate triplets. Each coordinate triplet can be defined either using an atom selection or three real numbers. The first vector is defined by coordinates P1 and P2, and the second vector, by coordinates P3 and P4. In addition, each vector can be defined in a different reference frame, specified optionally as `FR1<int>` and `FR2<int>`; the default frame is zero for each vector, which corresponds to the absolute coordinate frame. The angle and its derivatives are computed after the vectors are implicitly rotated into the same reference frame. The QUATERNION CV definition requires the specification of two coordinate frames FR1<int> and FR2<int>; the quaternion that is defined corresponds to the rotation of frame 2 into frame 1. Frames must be different (otherwise the quaternion is the identity rotation). -------------------------------------------------------------------------- [REPA] Call reparametriation (Repa) module. The supported reparametrization methods are the same as those described for the finite-temperature string module (FTSM) and the zero-temperature string (SM0K). -------------------------------------------------------------------------- [LIST] List currently defined collective variables. -------------------------------------------------------------------------- [CLEAr] Remove all collective variables, frames and quaternions. -------------------------------------------------------------------------- [TEST [GRAD {STEP <real>}|PARA|MINV] Test gradients of the string restraint potentials by finite difference if GRAD is specified; test parallel force computation when PARA is specified; test the computation of M-tensor inverse by LU decomposition, and by direct multi-diagonal matrix inversion. STEP sets the finite difference interval to use with GRAD. Parallel force computation test will work only when more than one CPU is assigned to a string image. The MINV test is provided for cross-validation. -------------------------------------------------------------------------- [FILL {MAIN|REF|OLD} {COMP} ] Calculate values of the CV from instantaneous coordinates and put into string column specified by one of "MAIN", "REF" or "OLD". This is how an initial string is defined in the beginning of a set of optimizations. Omission of column will result in the main CV coordinates being populated. Reference coordinates for statistics (see e.g. STAT RMSD) are assumed to be in the reference column (REF). The OLD column is used internally during string evolution, and contains the previous CV values (i.e. those prior to the most recent instance of 'EVOL'); the old values are used to linearly adjust CV coordinates from the old values to new values to avoid numerical instabilities. The user can also populate the OLD column manually, as can be useful in performing steered molecular dynamics (SMD). In this case, 'SMD' or 'STEERED' must be included in the 'STRING COLVAR DYNA' command. COMP indicates that CHARMM comparison coordinates are to be used for computing CV. -------------------------------------------------------------------------- [PRINt {ALL} {UNIT <int>} {NAME <character*>} {COL}] Write current CV coordinates to file from the column specified by COL. `ALL` indicates that the root node of each string replica is to write a separate file containing only the CV values of the corresponding replica. If ALL is omitted, one file will be created with the number of lines equal to the number of defined CV, each line containing M values, corresponding to M CVs. -------------------------------------------------------------------------- [READ {ALL | SCOL<int> TCOL<int>} {UNIT <int>} {NAME <character*>} {COL}] Read CV coordinates from file into the column specified by COL. `ALL` indicates that the root node of each string replica is to read a separate file containing only the CV values of the corresponding replica. If ALL is omitted, all CV coordinates are assumed to be contained in one with the number of lines equal to the number of defined CV, each line containing M values, corresponding to M CVs. If SCOL<i> is specified, each replica will read the CV coordinates of the i-th replica (i-th column) provided in the file (rather than it's corresponding column); TCOL must be provided along with SCOL, and indicates the number of columns (string images) in the file. -------------------------------------------------------------------------- [SWAP 2x[MAIN|REF|OLD]] Swap string coordinates in the specified columns. -------------------------------------------------------------------------- [COPY 2x[MAIN|REF|OLD]] Copy string coordinates from the first column to the second column. -------------------------------------------------------------------------- [SET [{IND<int> | ALL }] <parameter>] Will set parameter value for CV with index IND<i>, or for all CV if ALL is specified. <parameter> can be one or more of the following: [FORCe<real>] force constant for harmonic restraint potential. [GAMMa<real>] friction constant for CV evolution (see EVOL in DYNAmics below). [WEIGht<real>] CV weight for computing (string) distances. [ZVAL<real> [REP<int>] [{COL<int>}]] value for the CV variable, set for the specified replica in the specified column. -------------------------------------------------------------------------- [VOROnoi VCUT <real> {REP <int>} ] Set the maximum allowed distance from the string in the plane perpendicular to the string tangent that the instantaneous MD replica is allowed to travel in Voronoi simulations. This allows one to limit the transition tube width in Voronoi calculations. REP<i> will set the Voronoi cutoff distance on replica <i> only. [VOROnoi VMAP CALC] Compute the Voronoi map, which lists the Voronoi cell in which each instantaneous string coordinate set resides. In the present implementation, the functionality is used to ascertain that each instantaneous replica corresponds to the correct string image; i.e. MD replica on node (i) is closest to the string image (i), in which case the Voronoi map will consist ot two identical rows with entries 1 ... M, where M is the number of string images. Instantaneous coordinates and string coordinates must be defined prior to invoking this command. [VOROnoi VMAP PRIN [{UNIT <int>}] [NAME <character*>] ] Print Voronoi map to file with provided NAME and (optionally) provided unit number. [VOROnoi VMAP READ [{UNIT <int>}] [NAME <character*>] ] Read Voronoi map from file with provided NAME and (optionally) provided unit number. [VOROnoi VMAP CLEA] Clear the Voronoi map if it has been READ or CALCulated. This will force the Voronoi map to be recomputed at the beginning of dynamics. This option is important due to the present issue warnings if an MD replica is found outside its corresponding Voronoi cell; however, each replica is technically allowed to leave its cell for one iteration, at which point its (half-kick) momenta (the Verlet integrator is assumed) are reversed to return it to its home cell. In particular, is it possible that a replica will be outside of its cell at the timestep at which the restart file is written, in which case, the Voronoi map will be flagged as incorrect unless it is CLEARed prior to dynamics. [VOROnoi PRINt[{UNIT <int>}] [NAME <character*>] ] Print matrix C_ij with the number of Voronoi crossing attempts from cell (i) to cell (j) to file with provided NAME and (optionally) provided unit number. This file can be quickly postprocessed for a fast calculation of free energies of the tessellation. [VOROnoi READ [{UNIT <int>}] [NAME <character*>] ] Read matrix C_ij with the number of Voronoi crossing attempts from cell (i) to cell (j) to file with provided NAME and (optionally) provided unit number. This command is useful for restarting a Voronoi calculation with previously accumulated crossing statistics. -------------------------------------------------------------------------- [FRAMe <char*>] The FRAMe commands concern the manipulation of coordinate reference frames that can be used to compute rigid-body invariant positions, angles between vectors in different reference frames, and orientation quaternions. Addition of new frames is described under the ADD command above. [CLEAr] Remove all reference frames. [RESEt] Force any code that depends on a reference frame to recompute all reference frame axes. E.g., evaluating of a position CV expressed in a non-absolute frame of reference. [FILL [{COMP}] [{ALIGn}]] Compute reference frame axes from instantaneous coordinates according to the atom selection specified in ADD FRAME command. Because frame vectors are defined as eigenvectors of moment of inertia tensors of a set of atoms, frame vectors are only unique up to a sign. Considering only right-handed coordinate frames, there are three possible definitions of a frame. To ensure that different replicas define the frame vectors consistently, the option ALIGn is provided. When ALIGn is present, a reference coordinate set is assumed to be present in the comparison set, and the frame vectors computed from the main coordinate set are consistent with those computed from the comparison set. If the comparison set has the same coordinates for all string replicas, this almost surely implies that the frames are also defined consistently between the different replicas. COMP swaps the roles of the main and the comparison sets. [PRINt[{ALL}] [{UNIT<int>}] [{NAME<character*>}]] Write frame vectors to file with specified name and optional unit number. If ALL is specified, the root node of each string group/image will write a separate file with the frame vectors corresponding to the image. If ALL is omitted, frame vectors will be written for all replicas by the root node for the entire string. Each line in the files consists of nine entries corresponding to the three frame vectors (written out columnwise). The frame axes corresponding to different replicas are written on consecutive lines. If more than one frame is defined, the first frame is written out as a block of M lines (where M is the number of replicas); for N frames, there will be N such consecutive blocks. [READ [{ALL}] [{UNIT<int>}] [{NAME<character*>}]] Read frame vectors to file with specified name and optional unit number. File format and `ALL` option is as described for PRINt above. [LIST] List all defined frames and their constituent atoms. [ALIGn [{RMSD | VORO}]] Align frames along the string, making sure that the frames are defined consistently on different replicas. RMSD will cause the code to find the frame vectors such that the (position) CV values computed from the coordinates are closest to the string (in the notation of ref. [1] ||\theta(x) - z || is minimum). if VORO is specified, the code will select frame vectors such that the (position) CV values computed from the coordinates are closest to the string in the sense of the Voronoi metric, i.e. that each instantaneous coordinate set is within its corresponding Voronoi cell. The ALIGn oiptions should not normally be needed if the frame vectors are written out and read in during successive simulation, because in that case there is no ambiguity in their definition. -------------------------------------------------------------------------- [PARA [QUAT [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Indicate whether to enable parallel computation of quaternions and quaternion derivatives; default : off [PARA [FRAM [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Indicate whether to enable parallel computationl of reference frames and their derivatives; default : off [PARA [COLV [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Indicate whether to enable parallel computation of collective variables and their derivatives; default : off [PARA [MMAT [YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Indicate whether to enable parallel computation of the metric tensor; default : on [PARA [VORO[YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Indicate whether to enable parallel computation of Voronoi distances; default : on -------------------------------------------------------------------------- [MINV [LU|lu|DIAG|diag]] Specify method for inverting the metric tensor. LU specifies LU decomposition; DIAG specifies a multidiagonal direct inversion. DIAG is expected to be faster for matrices with nonzero entries clustered near the diagonal. The structure of the metric tensor depends on the definition of the CV. -------------------------------------------------------------------------- [MMAT [CALC]] Calculate the instantaneous metric tensor from coordinates. [MMAT PRINt[{INV}] [{UNIT<int>}] [NAME<character*>]] Write metric tensor to file with specified name and (optional) unit number. INV indicates that the inverse of the metric tensor is to be written. [MMAT READ [{INV}] [{UNIT<int>}] [NAME<character*>]] Read metric tensor from file with specified name and (optional) unit number. INV indicates that the inverse of the metric tensor is to be read. [MMAT [FAST[YES|ON|TRUE|T|yes|on|true|t|NO|OFF|FALSE|F|no|off|false|f]]] Specify whether to use a sparse matrix algorithm for computing the metric tensor. The FAST routine may not be faster if the metric tensor is not sparse (whether this is true depends on the definition of the CV). -------------------------------------------------------------------------- [WEIGht [CALC]] Compute CV weights from metric tensor. The weights are necessary do define distance in the space of CV. Alternatively, weights can be specified at the time of CV definition or read from file as described below. Whatever the method of choice, the weights should typically be held constant throughout a given SMCV simulation (e.g. they are equivalent to atomic masses if the CV are Cartesian coordinates). [WEIGht PRINt[{UNIT<int>}] [NAME<character*>]] Write CV weights to file with provided name and optional unit number. The output will consist of a single line with a number of floating point entries equal to the number of CV. [WEIGht READ [{UNIT<int>}] [NAME<character*>]] Read CV weights from file with provided name and optional unit number. -------------------------------------------------------------------------- [INTE [{INTERPCV [METH [LINE|SPLI|BSPL|LIN2]]}] - [NIN<int>] [CVIN<character*>] - [NOUT<int>][CVOUT<character*>] - [{COOR [CRIN<character*>][CROUT<character*>] - [{PDB [{RESI}] | FILE | UNFO | CARD | FORM}]]] [INTE] Interpolate a string onto a new string with different resolution. The routine can be used to interpolate CVs using the INTERPCV option, and also to create coordinate files that are optimal for the given CV values (starting from supplied coordinate files). [INTERPCV METH [LINE|SPLI|BSPL|LIN2] ] interpolation method to use. The valid methods are as described in the zero-temperature string documentation; default : LINE [NIN] Number of images in the existing string. [CVIN <character*>] Name of text file which contains the existing CV values (format given in PRINt command above). [NOUT] Number of coordinate sets in the interpolated string. [CVOUT <character*>] Name of text file to which interpolated CV values will be written. If INTERPCV is not set, the file specified by CVOUT is assumed to exist, containing the CV values for which coordinate files are to be generated (see CRIN/CROUT below). [{COOR <character*>}] When provided, coordinate files will be written in correspondence to the CV values in file CVOUT<char*>, using the following parameters. [CRIN <character*>] Name of text file that contains the filenames in which existing coordinates are stored, with one file name per line. [CROUT <character*> ] Name of text file that contains file names to which coordinate sets corresponding to the CVs in CVOUT<char*> are to be written. The following options are passed to the CHARMM coordinate input/output routines: [PDB [{RESI}] ] Indicates that input/output coordinate files are in PDB format. RESI will pass the 'RESID' option to the coordinate reader. [FILE | UNFO ] Indicates that input/output coordinate files are UNFORMATTED. [CARD | FORM}] Indicates that input/output coordinate files are FORMATTED. ------------------------------------------------------------------------- [DYNAmics] The STRIng COLV DYNAmics command will parse options relevant for the SMCV and subsequently call the regular dynamics integrator. Currently, SMCV support is limited to the default Verlet leapfrog integrator (dynamc module). The following, all optional, keywords are supported. [{STAF <int>}] Number of steps to skip between successive outputs of statistics. STAT options must be set prior to dynamics. [{NOPRint}] Reduce output from FTSM during dynamics. [{VORO}] For each FTSM replica, dynamics will be constrained to the corresponding Voronoi cell by reversing momenta at cell boundaries. This allows computing the free energy of the Voronoi tessellation and the MFPT between different Voronoi cells. [{EVOL [EVOF <int>] [{EVOS <int>}] [{EXPO {MEMO <real>}|AVER [NAVE <int>] }]}] Specify string evolution options. EVOL turns on string evolution during dynamics. EVOF <int> specifies the frequency at which the evolution is to be performed. A reasonable value for EVOF is 10, implying, essentially that every 10th structure during dynamics will be added to a running average. EVOS <int> specifies the number of dynamics steps after each update (if UPDA is set) during which evolution is temporarily turned off. This can be used to set an equilibration period after each update. A reasonable choice is to use the value for REEQ (see below) which specifies the number of dynamics steps over which the restraint potentials are adjusted to new string coordinates after an update. EVST <real> specifies the evolution step for the string images. By default, string evolution follows steepest descent dynamics on the free energy landscape of the collective variables (refs. 1-2), i.e., cv_i(n+1) = cv_i(n) - dt / gamma_i * M_ij * dF(cv_k; k\in[1,M])/dcv_j, where n is the nth iteration, dt is specified by EVST, gamma is the friction coefficient specified in the definition of each CV, M is the metric tensor, and F is the (M-dimensional) free energy as a function of the M collective variables. In the above equation, the quantity (dt/gamma) must units of squared time, where the unit of time is determined implicitly by the units of the MD engine; e.g. in CHARMM, the AKMA time unit is used. Reasonable values for dt (EVST) and gamma depend on the collective variable and the evolution frequency, but O(-3) and O(1), respectively, are a good place to start. EXPO {MEMO} specifies that evolution will use the exponential averaging, i.e., CV(N+1) = <rMEMO> *CV(N) + (1-<rMEMO>) * X_INST, where CV(N) are the CV coordinates at the Nth evolution call, X_INST are the instantaneous CV values computed from MD coordinates, and MEMO is a memory parameter in the range [0, 1]. Reasonable values for MEMO are 0.999 - 0.9999. AVER {NAVE <int>} Evolution will use simple averaging, i.e., PHI(N)=AVERAGE_{i=1}^{N}(X_i). To dampen oscillations at the beginning of a simulation, NAVE can be used to set an artificial number of samples in the average. That is, the average is updated using PHI(N+1) = ( NAVE * PHI(N) + X_INST ) / (NAVE++) ; therefore a high value for NAVE will decrease | PHI(N+1)-PHI(N) |. [{REPA [REPF <int>]}] indicates that reparametrization is to be performed after each block of REPF<int> iterations. REPA options must be set prior to dynamics. RSTR [REEQ <int>] {SMD|STEER} RSTR specifies that SMCV restraints are to be turned on during dynamics. REEQ specifies the number of simulation steps during which the restraint potentials will switch to the EVOLved CV coordinates from previous CV coordinates (during each evolution, if EVOL EVOF<i> is set). This option prevents abrupt changes in the restraints which can lead to instability. Values in the range 10 -- 500 are reasonable, depending on the system. When SMD or, equivalently, STEER are included, the OLD CV coordinate set (which can be FILLed manually, as described above) is not initialized prior to launching dynamics. This can be used to perform steered MD, whereby the CV values to which the MD replicas are restrained are computed by interpolating linearly between the OLD and MAIN CV sets. The steering is performed over the number of iterations set by REEQ. REX [REXF <int>] [REXT <real>] Instructs the code to attempt to swap MD coordinates between adjacent replicas at frequency REXF <i> using temperature REXT <r> in the Metropolis criterion. This option can significantly improve the convergence of free energy profiles but requires that the force constants in the restraining potentials be sufficiently low to facilitate frequent exchanges (a success rate of at least 10% is recommended). REXT can be increased beyond the simulation temperature to facilitate exchanges, but this will cause the thermodynamic ensemble to deviate from the Boltzmann distribution. ------------------------------------------------------------------------ [STATistics] Call statistics module. Specify statistics options when arguments are present. When called without arguments, an instance of statistics output will be written out. The SMCV STAT module has similar syntax to the zero-temperature (SM0K) and finite-temperature (FTSM) STAT modules. The specification of keywords RMSD RMSA DELS ARCL CURV FREE FORC VORO VMAP VLOG REXM REXL is as described for the SM0K and FTSM STAT routines. Keywords MMAT and COLV are described below. [STAT MMAT [MNAM <character*>]] Indicate that the metric tensor is to be output at each statistics call. The output file with the specified name will contain a column of N MxM matrices where N is the number of string replicas and M is the number of CV. [STAT COLV [CNAM <character*> {CVAP}]] Indicate that CV values are to be output at each statistics iteration to the fiel with the specified name, appending to an existing file (e.g. from previous calculations) is CVAP is set. Each line of the file corresponds to a single CV and contains N values for the N replicas. Successive lines correspond to different CV. ------------------------------------------------------------------------ Test cases. c40test/smcv.inp c40test/smcv-pos.inp c40test/smcv-voro.inp ------------------------------------------------------------------------ SMCV References: [1] L. Maragliano, A. Fischer, E. Vanden-Eijnden, and G. Ciccotti, J. Chem.Phys. 125, 024106 (2006). [2] L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006). [3] V. Ovchinnikov, M. Karplus, and E. Vanden-Eijnden, J. Chem. Phys. 134, 085103 (2011). If you find the SMCV code useful in preparing a manuscript, please cite the implementation reference 3 above.
Top Finite-Temperature String Method (FTSM). ----------------------------------------------------------------------------------------- Command Syntax. STRIng FTSM [INITialize] | [DONE] | [ADD [ORIEnt|RMSD] <atom selection>] | [LIST] | [CLEAr [ORIEnt|RMSD]] | [REPArametrize [{DEFI <real>}] - [{ITER <int>}] - [{LINE | SPLI | BSPL | DST [{WNCT <real>}] | LIN2}] | [ORIE <atom selection> {MASS}] | [TEST [GRAD {STEP <real>}|PARA] | [FILL {MAIN|REF} {COMP} ] | [COMM [ALLGather|HYPERcube|GATHERBroadcast]] | [PRINt [COR|DCD] {UNIT <int>} {NAME <character*>} {COL}] | [READ [COR|DCD] {UNIT <int>} {NAME <character*>} {COL}] | [SWAP 2x[MAIN|REF]] | [COPY 2x[MAIN|REF]] | [SET [ORIE <atom_selection>] | - [RMSD <atom_selection>] | - [KPAR <real> {ANGStrom}] | - [KPRP <real> {ANGStrom}] | - [KRMS <real> {ANGStrom}] | - [DPRP <real> {ANGStrom} {REP <int>} ] - [DRMS <real> {REP <int>} {UPPE} ] - [MASS {yes|on|true|t|YES|ON|TRUE|T|no|off|false|f|NO|OFF|FALSE|F}] - | [PROJ {yes|on|true|t|YES|ON|TRUE|T|no|off|false|f|NO|OFF|FALSE|F}] - | [VCUT <real> {REP <int>} ] ] | [VOROnoi [VMAP [CALC] | - [PRIN [{UNIT <int>}] [NAME <character*>] ] - [READ [{UNIT <int>}] [NAME <character*>] ] - [CLEA] ] | - [READ [UNIT <int>] [NAME <character*>] ] - [PRINt[UNIT <int>] [NAME <character*>] ] ] | [DYNAmics[{NOPR}] - [{VORO}] - [{UPDA [UPDF <int>] {REPA} ] - [{STAT [STAF]}] - [{EVOL [EVOF <int>] - [{EVOS <int>}] - [{EXPO {MEMO <real>} | AVER [NAVE <int>] {MAXAVE <int>}}]}]- [{RSTR [REEQ <int>]}] - [{REX [REXF <int>] [REXT <real>]}] - <regular dynamics options>] | [STATistics [{COUNt <int>]} - [{RMSD [RNAM <character*>] [{RAPP}]} ] - [{ARCL [ANAM <character*>] [{AAPP}]} ] - [{CURV [CVNM <character*>] [{CAPP}]} ] - [{DIST [DNAM <character*>] [{DAPP}]} ] - [{FREE [FENM <character*>] [{FAPP}] [{NOCV}] }] - [{FORC [FCNM <character*>] [{FCAP}] }] - [{CENT [CNAM <character*>] [{CEAP}] }] [{REXM [RXNM <character*>] [{RXOL <character*>}]}] - [{REXL [{RXNM <character*>}] [{ROFF <int>}] [{RXAP}]}] - [{VORO [VNAM <character*>]} ] - [{VMAP [{VNAM <character*>}]} ] - [{VLOG [{VNAM <character*>}]} [{VOFF <int>}] [{VLAP}]}] Command Description. -------------------------------------------------------------------------- [INIT] Initialize the finite-temperature string module. Substitution parameters ?MESTRING and ?NSTRING will be set to correspond to the string rank of each CPU group, and to the total number of string replicas, respectively. -------------------------------------------------------------------------- [DONE] Finalize the finite-temperature string module -------------------------------------------------------------------------- [REPA] Call reparametrization (Repa) module. The supported reparametrization methods are the same as those described for the zero-temperature string method (SM0K), except that keywords ORIE, MASS and MOVE are not supported (similar functionality is configured by other means described below) -------------------------------------------------------------------------- [ADD [ORIEnt|RMSD] <atom selection>] Add a center-of-mass (COM) position vector computed from the atom selection to the reaction set (i.e. coordinates that define the string). ORIEnt will add to the orientation set (groups which will define the best-fit superposition); RMSD will add to the forcing set (groups between which the distances are computed after alignment. NOTE: an empty orientation set implies that best-fit alignments will not be performed. COM groups should be nonoverlapping, otherwise the metric tensor will not, generally, be close to the corresponding inverse mass tensor (see ref. [1] for a discussion of the metric tensor), which will probably lead to inaccurate free energies; however this condition is not checked in the code. -------------------------------------------------------------------------- [LIST] List currently defined orientation and forcing sets. -------------------------------------------------------------------------- [CLEAr [ORIEnt|RMSD]] Clear orientation or forcing set. -------------------------------------------------------------------------- [ORIE <atom selection> {MASS}] Perform best-fit superposition of the instantaneous coordinates, such that the root-mean-squared distance (RMSD) between instantaneous coordinates of string image (i>0) and (i-1) are minimum. MASS specifies that mass-weighting is to be used. -------------------------------------------------------------------------- [TEST [GRAD {STEP <real>}|PARA] Test gradients of the string restraint potentials by finite difference if GRAD is specified; test parallel force computation when PARA is specified. STEP sets the finite difference interval to use with GRAD. Parallel force computation test will work only when more than one CPU is assigned to a string image. -------------------------------------------------------------------------- [FILL {MAIN|REF} {COMP} ] Calculate values of the COM positions from instantaneous coordinates and put into string column specified by MAIN or REF. This is how an initial string is defined in the beginning of a set of optimizations. Omission of column will result in the main column being populated. Reference coordinates for statistics (see e.g. STAT RMSD) are assumed to be in column REF. COMP specifies that CHARMM comparison coordinates are to be used for computing COM positions. -------------------------------------------------------------------------- [COMM [ALLGather|HYPERcube|GATHERBroadcast]] Specify parallel communication scheme for gathering gradients computed in parallel. ALLG uses MPI_allgather ; HYPER uses an internal hypercube subroutine implemented via MPI_SendRecv ; GATHERB uses MPI_Gather to local root followed by MPI_Bcast. This choice should not normally have a large impact on performance. The default method is hypercube. -------------------------------------------------------------------------- [PRINt [COR|DCD] {UNIT <int>} {NAME <character*>} {COL}] Write current string coordinates to file from the column specified by COL. COR and DCD specify CHARMM coordinate and trajectory format, respectively. If COM groups contain more than one atom, the coordinate of the COM group will be output as the coordinate of the first atom in the COM group. It is possible to define multiple COM groups with the same first atom. However, this possibility is not considered at present (see short discussion above following ADD). -------------------------------------------------------------------------- [READ [COR|DCD] {UNIT <int>} {NAME <character*>} {COL}] Read string coordinates from file into the column specified by COL. COR and DCD specify CHARMM coordinate and trajectory format, respectively. If COM groups contain more than one atom, the coordinate of the COM group will be found as the coordinate of the first atom in the COM group. -------------------------------------------------------------------------- [SWAP 2x[MAIN|REF]] Swap string coordinates in the specified columns. -------------------------------------------------------------------------- [COPY 2x[MAIN|REF]] Copy string coordinates from the first column to the second column. -------------------------------------------------------------------------- [SET ORIE <atom_selection>] In addition to the ADD command for (iteratively) defining atom groups (see above), this command is a shortcut for defining the entire ORIEntation set, albeit with one atom per group, i.e. the number of groups is the number of atoms in the atom selection, and the COM position of the group is the position of the corresponding (single) atom in the group. [SET RMSD <atom_selection>] Same as above, but defines the forcing group [SET KPAR <real> {ANGStrom}] Set force constant for plane restraint (i.e. force acting parallel to the path; see ref. [1] for definition of restraint potentials). ANGStrom indicates that the distance unit in the force constant are Angstrom. If unspecified, the distance unit is equal to twice the distance between adjacent string replicas. [SET KPRP <real> {ANGStrom}] Set force constant for in-plane restraint (i.e. force acting perpendicular to the path for limiting). This option is used for limiting the transition tube width. ANGStrom option as above. [SET KRMS <real> {ANGStrom}] Set force constant for a regular RMSD restraint, whereby the MD replica is restrained to the string image via U = krms/2 || RMSD (x , x0) - RMSD0 ||^2 (where x includes an implicit best-fit rotation). This potential corresponds approximately to the behavior in the replica path functionality. Although the free energies cannot be computed, the RMSD potential is useful for path equilibration. ANGStrom option as above. [SET DPRP <real> {ANGStrom} {REP <int>} ] Set the reference distance for the in-plane restraint perpendicular to the path; used for limiting the transition tube width, as the MD replica will be restrained to a vicinity of the path. REP allows the use of different values for different replicas, although this flexibility is not often used. ANGStrom option as above. [SET DRMS <real> {REP <int>} {UPPE} ] Set the reference RMSD0 for the RMSD restraint above; if omitted, the default value of zero is used. If UPPE is specified, the restraint is only active when the computed RMSD(x) exceeds the reference value. I.e. this corresponds to the one-sided RMSD potential U = krms/2 max(0,|| RMSD (x , x0) - RMSD0||)^2 [SET MASS {yes|on|true|t|YES|ON|TRUE|T|no|off|false|f|NO|OFF|FALSE|F}] When set to on, all distances will be calculated after mass-weighting the string coordinates; best-fit superposition will also use mass weights. When off, uniform weights will be used. [SET PROJection {yes|on|true|t|YES|ON|TRUE|T|no|off|false|f|NO|OFF|FALSE|F}] When set to on, coordinate vectors pointing from the string image to the instantaneous MD replica are projected onto the path tangent. This allows the imposition of restraints to the distance to the hyperplane normal to the string. This option should normally be on for FTSM to allow the computation of free energies. The force constants corresponding to the parallel and perpendicular force are KPAR and KPRP defined above. When set to off, the MD replica is restrained to the corresponding image using RMSD restraints (see KRMS above). [SET VCUT <real> {REP <int>} ] Set the maximum allowed distance from the string in the plane perpendicular to the string tangent that the instantaneous MD replica is allowed to travel in Voronoi simulations. This allows one to limit the transition tube width in Voronoi calculations, and is the constrained counterpart of DPRP/KPRP above. REP<i> will set the Voronoi cutoff distance on replica <i> only. ------------------------------------------------------------------------ [VOROnoi VMAP CALC] Compute the Voronoi map, which lists the Voronoi cell in which each instantaneous string coordinate set resides. In the present implementation, the functionality is used to ascertain that each instantaneous replica corresponds to the correct string image; i.e. MD replica on node (i) is closest to the string image (i), in which case the Voronoi map will consist ot two identical rows with entries 1 ... N, where N is the number of string images. Instantaneous coordinates and string coordinates must be defined prior to invoking this command. [VOROnoi VMAP PRIN [{UNIT <int>}] [NAME <character*>] ] Print Voronoi map to file with provided NAME and (optionally) provided unit number. [VOROnoi VMAP READ [{UNIT <int>}] [NAME <character*>] ] Read Voronoi map from file with provided NAME and (optionally) provided unit number. [VOROnoi VMAP CLEA] Clear the Voronoi map if it has been READ or CALCulated. This will force the Voronoi map to be recomputed at the beginning of dynamics. This option is important due to the present issue warnings if an MD replica is found outside its corresponding Voronoi cell; however, each replica is technically allowed to leave its cell for one iteration, at which point its (half-kick) momenta (the Verlet integrator is assumed) are reversed to return it to its home cell. In particular, is it possible that a replica will be outside of its cell at the timestep at which the restart file is written, in which case, the Voronoi map will be flagged as incorrect unless it is CLEARed prior to dynamics. [VOROnoi PRINt[{UNIT <int>}] [NAME <character*>] ] Print matrix C_ij with the number of Voronoi crossing attempts from cell (i) to cell (j) to file with provided NAME and (optionally) provided unit number. This file can be quickly postprocessed for a fast calculation of free energies of the tessellation. [VOROnoi READ [{UNIT <int>}] [NAME <character*>] ] Read matrix C_ij with the number of Voronoi crossing attempts from cell (i) to cell (j) to file with provided NAME and (optionally) provided unit number. This command is useful for restarting a Voronoi calculation with previously accumulated crossing statistics. ------------------------------------------------------------------------ [DYNAmics] The STRIng FTSM DYNAmics command will parse options relevant for the FTSM and subsequently call the regular dynamics integrator. Currently, FTSM support is limited to the default Verlet leapfrog integrator (dynamc module). The following, all optional, keywords are supported. [{STAF <int>}] Number of steps to skip between successive outputs of statistics. STAT options must be set prior to dynamics. [{NOPRint}] Reduce output from FTSM during dynamics. [{VORO}] For each FTSM replica, dynamics will be constrained to the corresponding Voronoi cell by reversing momenta at cell boundaries. This allows computing the free energy of the Voronoi tessellation and the MFPT between different Voronoi cells. [{UPDA [UPDF <int>] {REPA}] When UPDA is present string coordinates will be updated from the running averages computed by option EVOL. UPDF specifies the frequency of updating (i.e. number of the intervening dynamics steps). REPA indicates that a raperametrization is to be performed after each update. REPA options must be set prior to dynamics. [{EVOL [EVOF <int>] [{EVOS <int>}] [{EXPO {MEMO <real>}|AVER [NAVE <int>] {MAXAVE <int>}}]}] Specify string evolution options. EVOL turns on string evolution during dynamics. EVOF <int> specifies the frequency at which the evolution is to be performed. A reasonable value for EVOF is 10, implying, essentially that every 10th structure during dynamics will be added to a running average. EVOS <int> specifies the number of dynamics steps after each update (if UPDA is set) during which evolution is temporarily turned off. This can be used to set an equilibration period after each update. A reasonable choice is to use the value for REEQ (see below) which specifies the number of dynamics steps over which the restraint potentials are adjusted to new string coordinates after an update. EXPO {MEMO} specifies that evolution will use the exponential averaging, i.e., PHI(N+1) = <rMEMO> *PHI(N) + (1-<rMEMO>) * X_INST, where PHI(N) are the string coordinates at the Nth evolution call, X_INST are the instantaneous MD coordinates (after best-fit superposition, if the orientation set is non-empty), and MEMO is a memory parameter in the range [0, 1]. Reasonable values for MEMO are 0.999 - 0.9999. EXPO AVER {NAVE <int>} {MAXAVE <int>} Evolution will use simple averaging, i.e., PHI(N)=AVERAGE_{i=1}^{N}(X_i). To dampen oscillations at the beginning of a simulation, NAVE can be used to set an artificial number of samples in the average. That is, the average is updated using PHI(N+1) = ( NAVE * PHI(N) + X_INST ) / (NAVE++) ; therefore a high value for NAVE will decrease | PHI(N+1)-PHI(N) |. Specifying MAXAVE <i> will prevent NAVE from exceeding <i>; this effectively transforms this evolution method to EXPOnential evolution described above. MAXAVE is zero by default. RSTR [REEQ <int>] RSTR specifies that FTSM restraints are to be turned on during dynamics. The type of restraint is determined by PROJection [on|off] (see above). REEQ specifies the number of simulation steps during which the restraint potentials will switch to using UPDAted string coordinates from previous string coordinates (during each update, if updates are on). This option prevents abrupt changes in the restraints which can lead to instability. Values in the range 10 -- 500 are reasonable, depending on the system. REX [REXF <int>] [REXT <real>] Instructs the code to attempt to swap MD coordinates between adjacent replicas at frequency REXF <i> using temperature REXT <r> in the Metropolis criterion. This option can significantly improve the convergence of free energy profiles but requires that the force constants in the restraining potentials be sufficiently low to facilitate frequent exchanges (a success rate of at least 10% is recommended). REXT can be increased beyond the simulation temperature to facilitate exchanges, but this will cause the thermodynamic ensemble to deviate from the Boltzmann distribution. ------------------------------------------------------------------------ [STATistics] Call statistics module. Specify statistics options when arguments are present. When called without arguments, an instance of statistics output will be written out. The FTSM STAT module has similar syntax to the zero-temperature (SM0K) STAT module. [COUNt <int>] specifies iteration counted for statistics output. The default value is 0. This number is incremented after every call to statistics, and is printed in the output files corresponding to RMSD, ARCL, CURV, DIST, FREE and FORCE. [RMSD RNAM <character*> [RAPP]] sets up output of RMSD values between the replica coordinates at the current iteration and the coordinates present in the comparison set (usually the initial string). At each call to string statistics, a line will be added to the file with the name specified in [RNAM <character*>], with the columns in the file corresponding to different replicas. If RAPP is specified, output will be appended to the file. If RNAM is omitted, RMSD output is directed to the output stream. This subcommand is useful for gauging convergence of the string. [ARCL ANAM <character*> [AAPP]] sets up output of the distance between the adjacent replicas (such that their sum yields the string length) at the current iteration. At each call to string statistics, a line will be added to the file with the name specified with [ANAM <character*>]. If AAPP is specified, output will be appended to the file. [CURV CVNM <character*> [CAPP]] sets up output of the approximate string curvature at each string image (the curvature is approximate because it is computed under the assumption that the distance between adjacent string replicas is exactly equal). At each call to string statistics, a line will be added to the file with the name specified with [CVNM <character*>]. If CAPP is specified, output will be appended to the file. NOTE: arclength and curvature are computed by the reparametrization routine; therefore REPA must be active during dynamics to output ARCL and CURV statistics (otherwise the corresponding files will have zeros). If reparametrization is not desired, it can be used with 'REPA ITER 0', which will force the routine to be called but will not allow it to perform any reparametrization iterations. DIST [DNAM <character*>] [{DAPP}]] sets up output of the coordinate projections onto the string at each string image. At each call to string statistics, one or two lines will be added to the file with the name specified with [DNAM <character*>], depending in whether PROJ is on|off. If DAPP is specified, output will be appended to the file. The output can be understood from the following diagram r /| }perpendicular distance ---p---o---q--- \_____/ d -- projection onto q-p \_______/ D -- displacement vector q-p where r is the instantaneous coordinate assigned to string image o, p and q are the adjacent string images. If PROJ is on, the first line is the normalized projection of r-o onto q-p, i.e. (r-p) . (q-p) ------------- = d/||D|| \def delta ||q-p||^2 and the second line is the normalized distance to the string, i.e. || ( r-p - d(q-p)/||D||) || / ||D|| These are the variables that are restrained in FTSM calculations when PROJ is on. For example, for internal string points, delta is restrained to 0.5, and for the left and right endpoints, to 0 and 1, respectively. Note that ||D|| equals twice the string length increment ds. If PROJ is off, ||r-o|| will be printed. In the above discussion, best-fit rotations are implicit when the orientation set is nonempty. FREE [FENM <character*>] [{FAPP}] [{NOCV}] sets up output of the free energy (FE). At each call to string statistics, a line will be added to the file with the name specified with [FENM <character*>]. If FAPP is specified, output will be appended to the file. If NOCV is specified, curvature contributions will not be included in the FE output. FORC [FCNM <character*>] [{FCAP}] sets up output of the thermodynamic forces that can be integrated to obtain the free energy. If PROJ is on, at each call to string statistics, three lines will be added to the file with the name specified with [FENM <character*>] (only two lines are added if FREE energy output is requested as above and NOCV is specified). If FAPP is specified, output will be appended to the file. The first line contains the forces acting on the parallel projection variable (delta above). The ds metric is already included in the force, so the FE at replica (i) can be computed (to first order) simply as (sum force_{j=1}^i ), although second-order (or higher) integration schemes such as the trapezoid rule are preferable. The second line contains the curvature contribution to the force (-\xi} in ref [1]. It can be added to the projection force to account for the curvature, if desired. If NOCV is specified with FREE energy output (above) this line will be omitted. The third line lists the forces acting on the distance from the string used to limit transition tube width. If is provided for information ans is not useful for computing the FE. If PROJ is off, the force acting on the RMSD variable (see section in KRMS above). This force is provided for information and cannot be used to compute the free energy. CENT [CNAM <character*>] [{CEAP}] sets up output of the string images in trajectory (.dcd) format. At each statistics call N coordinate sets will be written consecutively to the specified file, where N is the number of string images. If COM groups contain more than one atom, the coordinate of the COM group will be output as the coordinate of the first atom in the COM group. It is possible to define multiple COM groups with the same first atom. However, this possibility is not considered at present (see short discussion above following ADD). REXM [RXNM <character*>] [{RXOL <character*>}] sets up output of the replica exchange map, which the current string image to which each MD replica is currently assigned (in the replica exchange functionality, the MD replica coordinates are allowed to `travel` between string images according to the Metropolis criterion to facilitate equilibration) RXNM is the name of the replica map file to which extension .map will be appended. RXOL allows the specification of an existing replica map (e.g. from a previous dynamics run) so that replica jump dynamics can be monitored continuously across restarts. REXL [{RXNM <character*>}] [{ROFF <int>}] [{RXAP}] sets up output of the replica exchange log, which contains a record of all the successful swap attempts. Each line lists the timestep at which the jump was performed and the two replicas whose coordinates have been swapped. ROFF is an offset added to the timestep counter to provide continuous time evolution from previous runs. RXAP indicates that an existing log file is to be appended to. VORO [VNAM <character*>] specify name of file to which to print the matrix of crossing attempts C_ij (see VOROnoi commands above). The extension '.dat' will be appended to the file name. VMAP [{VNAM <character*>}] specify name of file to which to print the Voronoi map (see VOROnoi commands above). The extension '.map' will be appended to the file name. VLOG [{VNAM <character*>}] [{VOFF <int>}] [{VLAP}] specifies the name of the binary Voronoi log file which contains the time record of all crossing attempts. This file is written in binary format because it can become very large in order to accumulate sufficient statistics to compute free energies and mean first passage times. The file is post- processed using the MATLAB/Octave script vcalc.m given in the support directory. Briefly, each record in the file contains a variable number of entries, each entry being composed of 5 integers: (1) the iteration (or timestep) number, (2) string replica ID, (3) ID of the Voronoi cell from which the crossing attempt is being made, (4) the ID of the Voronoi cell to which the crossing attempt is being made, (5) the ID of the cell in which the replica resides after crossing attempt (i.e. whether the replica was reflected or allowed to cross). The present version of the code does not allow crossing; every replica is reflected, and all crossing attempts `fail`; however, this need not be the case more generally, and this feature may be included in future releases. From the time record of crossing attempts, the free energy and reaction rates can be approximated according to refs. [3-4]. VOFF provides an optional offset for the timestep to preserve time-continuity of the crossing dynamics. VLAP indicates that an existing log file is to be appended to (however, the postprocessing script vcalc.m can accept multiple successive log files). The extension '.log' will be appended to the file name specified by VNAM. NOTE : The base file name specified by VNAM is shared by the the VORO output options above, and needs to be specified only once. ------------------------------------------------------------------------ Test cases. c40test/ftsm.inp c40test/ftsm-fix.inp c40test/ftsm-voro.inp ------------------------------------------------------------------------ References [1] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B. 66, 052301 (2002). [2] E. Vanden-Eijnden and M. Venturoli, J. Chem. Phys. 131, 044120 (2009). [3] V. Ovchinnikov and M. Karplus, J. Chem. Phys. 140, 175103 (2014). If you find the FTSM code useful in preparing a manuscript, please cite the implementation reference 3 above.
Top Test cases The following test cases apply the string methods to the c7eq/c7ax transition in the alanine dipeptide in vacuum. test/c40test/sm0k.inp Zero-temperature string test/c40test/smcv.inp String in phi/psi dihedral angles test/c40test/smcv-pos.inp String in Cartesian COM positions test/c40test/smcv-voro.inp String with tessellations test/c40test/ftsm.inp Finite-temperature string test/c40test/ftsm-fix.inp Finite-temperature string without best-fit alignment test/c40test/ftsm-voro.inp Finite-temperature string with tessellations c40test/confcons.inp Test of conformational consistency using a myosin 6 fragment Supporting files The following are matlab/octave files that can be used to postprocess and display output from string calculations. support/stringm/arc.m : plot string length support/stringm/rmsd.m : plot rmsd between evolving string and initial string support/stringm/smooth2.m : optional smoothing routine support/stringm/vcalc.m : compute and plot free energy and mean first passage time from voronoi logs support/stringm/voro_fe.m : compute and plot free energy using collision data from voronoi simulations support/stringm/ftsm/free.m : compute free energy from ftsm simulation data support/stringm/sm0k/ener.n : plot potential energies along zero-temperature string in sm0k data support/stringm/smcv/free.m : plot free energy output from SMCV simulations
Top Zero-temperature string method (SM0K) [1] W. E, W. Ren, and E. Vanden-Eijnden, J. Chem. Phys. 126, 164103 (2007). [2] V. Ovchinnikov, M. Karplus, and E. Vanden-Eijnden, J. Chem. Phys. 134, 085103 (2011). [3] I. Khavrutskii, K. Arora, and C. Brooks III, J. Phys. Chem. 125, 174108 (2006). If you find the SM0K code useful in preparing a manuscript, please cite the implementation reference 2 above. Finite-temperature string method (FTSM) [1] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B. 66, 052301 (2002). [2] E. Vanden-Eijnden and M. Venturoli, J. Chem. Phys. 131, 044120 (2009). [3] V. Ovchinnikov and M. Karplus, J. Chem. Phys. 140, 175103 (2014). If you find the FTSM code useful in preparing a manuscript, please cite the implementation reference 3 above. String method in collective variables (SMCV) [1] L. Maragliano, A. Fischer, E. Vanden-Eijnden, and G. Ciccotti, J. Chem.Phys. 125, 024106 (2006). [2] L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006). [3] V. Ovchinnikov, M. Karplus, and E. Vanden-Eijnden, J. Chem. Phys. 134, 085103 (2011). If you find the SMCV code useful in preparing a manuscript, please cite the implementation reference 3 above.