diff --git a/doc/colvars-code-refs.bib b/doc/colvars-code-refs.bib index 6aa885010..d4252bbd1 100644 --- a/doc/colvars-code-refs.bib +++ b/doc/colvars-code-refs.bib @@ -102,9 +102,9 @@ @article {Ebrahimi2022 % Moving frame of reference % Optimal rotation via flexible fitting % orientationAngle colvar component (derived from orientation) -% orientationProj colvar component (derived from orientation) -% spinAngle colvar component (derived from orientation) -% tilt colvar component (derived from orientation) +% orientationProj colvar component (derived from orientationAngle) +% spinAngle colvar component (derived from tilt) +% tilt colvar component (derived from orientationProj) % alpha colvar component % dihedralPC colvar component % cartesian colvar component @@ -161,9 +161,9 @@ @article{Fu2016 % polarTheta colvar component % polarPhi colvar component -% eulerPhi colvar component (derived from orientation) -% eulerTheta colvar component (derived from orientation) -% eulerPsi colvar component (derived from orientation) +% eulerPhi colvar component (derived from orientation_angle) +% eulerTheta colvar component (derived from orientation_angle) +% eulerPsi colvar component (derived from orientation_angle) @article{Fu2017, author = {Fu, Haohao and Cai, Wensheng and H\'enin, J\'er\^ome and Roux, Beno\^it and Chipot, Christophe}, title = {New Coarse Variables for the Accurate Determination of Standard Binding Free Energies}, diff --git a/doc/colvars-refman-main.tex b/doc/colvars-refman-main.tex index 709c897b1..d393b6cb7 100644 --- a/doc/colvars-refman-main.tex +++ b/doc/colvars-refman-main.tex @@ -97,49 +97,49 @@ The values of ``$c$''\cvnamebasedonly{ and ``alpha''} are also recorded throughout the simulation as a joint 2-dimensional histogram (see \ref{sec:colvarbias_histogram}). \begin{cvexampleinput} -colvar \{\\ -\-~~\# difference of two distances\\ -\-~~name d \\ -\-~~width 0.2 \# estimated fluctuation width \\ -\-~~distance \{\\ -\-~~~~componentCoeff 1.0\\ -\-~~~~group1 \{ atomNumbers 1 2 \}\\ -\-~~~~group2 \{ atomNumbers 3 4 5 \}\\ +colvar~\{\\ +\-~~\#~difference~of~two~distances\\ +\-~~name~d~\\ +\-~~width~0.2~~\#~estimated~fluctuation~width~\\ +\-~~distance~\{\\ +\-~~~~componentCoeff~~1.0\\ +\-~~~~group1~\{~atomNumbers~1~2~\}\\ +\-~~~~group2~\{~atomNumbers~3~4~5~\}\\ \-~~\}\\ -\-~~distance \{\\ -\-~~~~componentCoeff -1.0\\ -\-~~~~group1 \{ atomNumbers 7 \}\\ -\-~~~~group2 \{ atomNumbers 8 9 10 \}\\ +\-~~distance~\{\\ +\-~~~~componentCoeff~-1.0\\ +\-~~~~group1~\{~atomNumbers~7~\}\\ +\-~~~~group2~\{~atomNumbers~8~9~10~\}\\ \-~~\}\\ \}\\ \\ -colvar \{\\ -\-~~name c\\ -\-~~coordNum \{\\ -\-~~~~cutoff 6.0\\ -\-~~~~group1 \{ atomNumbersRange 1-10 \}\\ -\-~~~~group2 \{ atomNumbersRange 11-20 \}\\ -\-~~~~tolerance 1.0e-6\\ -\-~~~~pairListFrequency 1000\\ +colvar~\{\\ +\-~~name~c\\ +\-~~coordNum~\{\\ +\-~~~~cutoff~6.0\\ +\-~~~~group1~\{~atomNumbersRange~~1-10~\}\\ +\-~~~~group2~\{~atomNumbersRange~11-20~\}\\ +\-~~~~tolerance~1.0e-6\\ +\-~~~~pairListFrequency~1000\\ \-~~\}\\ \}\\ \cvnamebasedonly{% -colvar \{\\ -\-~~name alpha\\ -\-~~alpha \{\\ -\-~~~~psfSegID PROT\\ -\-~~~~residueRange 1-10\\ +colvar~\{\\ +\-~~name~alpha\\ +\-~~alpha~\{\\ +\-~~~~psfSegID~PROT\\ +\-~~~~residueRange~1-10\\ \-~~\}\\ \}} \\ -harmonic \{\\ -\-~~colvars d c\\ -\-~~centers 3.0 4.0\\ -\-~~forceConstant 5.0\\ +harmonic~\{\\ +\-~~colvars~d~c\\ +\-~~centers~3.0~4.0\\ +\-~~forceConstant~5.0\\ \}\\ \\ -histogram \{\\ -\-~~colvars c\cvnamebasedonly{ alpha}\\ +histogram~\{\\ +\-~~colvars~c\cvnamebasedonly{~alpha}\\ \} \end{cvexampleinput} @@ -313,7 +313,7 @@ Such a custom variable can either be computed on the fly (by calling an external tool or a VMD feature) or pre-computed and cached. For example, we can use VMD to compute the protein solvent-accessible surface area (SASA) on the fly: -\begin{cvexampleinput} +\begin{mdexampleinput} \-set~sasa\_sel~[atomselect~top~"protein~and~noh"]\\ \-\\ \-proc~calc\_user\_colvar~\{~x~\}~\{\\ @@ -321,7 +321,7 @@ \-~~\$::sasa\_sel~frame~[cv~frame]\\ \-~~return~[measure~sasa~1~\$::sasa\_sel]\\ \-\} -\end{cvexampleinput} +\end{mdexampleinput} The Tcl script above should be entered into VMD's text console, or sourced from a file. Then we can define a colvar that will be calculated by the scripted procedure above, by entering the following configuration as a new colvar in the Dashboad Editor Window: @@ -344,7 +344,7 @@ If the SASA computation is slow, it can be cached to make subsequent interaction responsive. To that effect, we source the following script to pre-compute the SASA for every frame and store it in VMD's User field for atom 0. -\begin{cvexampleinput} +\begin{mdexampleinput} \-set~atom0\_sel~[atomselect~top~"index~0"]\\ \-\\ \-for~\{set~f~0\}~\{\$f~<~[molinfo~top~get~numframes]\}~\{incr~f\}~\{\\ @@ -352,18 +352,18 @@ \-~~\$::atom0\_sel~frame~\$f\\ \-~~\$::atom0\_sel~set~user~[measure~sasa~1~\$::sasa\_sel]\\ \} -\end{cvexampleinput} +\end{mdexampleinput} Then we use the following Tcl procedure to compute a custom variable that is simply equal to the User field for atom 0 for the current frame: -\begin{cvexampleinput} +\begin{mdexampleinput} \-\\ \-proc~calc\_user\_colvar~\{~x~\}~\{\\ \-~~\#~Follow~current~Colvars~frame\\ \-~~\$::atom0\_sel~frame~[cv~frame]\\ \-~~return~[\$::atom0\_sel~get~user]\\ \-\} -\end{cvexampleinput} +\end{mdexampleinput} This version results in extremely fast updates. } % cvvmdonly dashboard @@ -851,7 +851,7 @@ \noindent{}where each argument is a string passed to the function or functions that are used to compute the variable, and are called colvar components, or CVCs (\ref{sec:cvc_list}). For example, a variable ``\texttt{DeltaZ}'' made of a single ``\texttt{distanceZ}'' component can be made periodic with a period equal to the unit cell dimension along the $Z$-axis: \cvscriptexampleinputcolvar{modifycvcs}{"period \$\{Lz\}"}{"DeltaZ"} -\noindent{}Please note that this option is currently limited to changing the values of the polynomial superposition parameters \refkey{componentCoeff}{sec:cvc_superp}, or of the \refkey{componentExp}{sec:cvc_superp} to update on the fly, of \refkey{period}{sec:cvc_periodic}, \refkey{wrapAround}{sec:cvc_periodic} or \refkey{forceNoPBC}{colvar|distance|forceNoPBC} options for components that support them. +\noindent{}Please note that this option is currently limited to changing the values of the polynomial superposition parameters \refkey{componentCoeff}{sec:cvc_superp}, or of the \refkey{componentExp}{sec:cvc_superp} to update on the fly, of \refkey{period}{sec:cvc_periodic}, \refkey{wrapAround}{sec:cvc_periodic} or \refkey{forceNoPBC}{colvar|cvc|forceNoPBC} options for components where it is relevant. \cvlammpsonly{Furthermore, because the call above uses an \emph{immediate} expansion of the variable ``\texttt{\$\{Lz\}}'', its value will not be updated during a constant-pressure simulation.} If the variable is computed using many components, it is possible to selectively turn some of them on or off: @@ -1358,7 +1358,7 @@ A collective variable is defined by the keyword \texttt{colvar} followed by its configuration options contained within curly braces: \begin{cvexampleinput} -\-colvar \{\\ +\-colvar~\{\\ \-~~name~xi\\ \-~~$<$other~options$>$\\ \-~~function\_name~\{\\ @@ -1476,6 +1476,37 @@ The very few keywords that are available for all types of components are listed in a separate section \ref{sec:cvc_common}. +\cvsubsec{Treatment of periodic boundary conditions}{sec:cvc_pbcs} + +In all colvar components described below, the following rules apply concerning periodic boundary +conditions (PBCs): +\begin{enumerate} +\item Distance vectors between two coordinates + $\mathbf{d}_{i,j} = \left(\mathbf{x}_1 - \mathbf{x}_2\right)$, are calculated following the + minimum-image convention by default, unless \refkey{forceNoPBC}{colvar|cvc|forceNoPBC} is enabled. + ($\mathbf{x}_1$ and $\mathbf{x}_2$ may be either individual atomic coordinates, or centers of mass + of two groups.) +\item For all other functions of individual atomic coordinates, + $f\left(\mathbf{x}_1, \mathbf{x}_2, \ldots\right)$, it is assumed that all atoms that are part of + the same group are in the same periodic unit cell (see \ref{sec:colvar_atom_groups_wrapping}). +\end{enumerate} + +\begin{itemize} +\item % + \labelkey{colvar|cvc|forceNoPBC}% + \keydef{forceNoPBC}{% + any component}{% + Use absolute rather than minimum-image distances?}{% + boolean}{% + \texttt{no}}{% + If this keyword is set to \texttt{yes}, PBCs will be ignored when calculating this component. + This is only useful in a limited number of special cases, e.g.\ to describe the distance between + remote points of a single macromolecule, which cannot be split across periodic cell boundaries, + and for which the minimum-image distance might give the wrong result because of a relatively + small periodic cell.} +\end{itemize} + + \cvsubsec{Distances}{sec:cvc_distances} @@ -1498,24 +1529,6 @@ \labelkey{colvar|distance|group2} \simkey{group2}{\texttt{distance}}{group1} -\item % - \labelkey{colvar|distance|forceNoPBC} - \keydef - {forceNoPBC}{% - \texttt{distance}}{% - Calculate absolute rather than minimum-image distance?}{% - boolean}{% - \texttt{no}}{% - By default, in calculations with periodic boundary conditions, the - \texttt{distance} component returns the distance according to the - minimum-image convention. If this parameter is set to \texttt{yes}, - PBC will be ignored and the distance between the coordinates as maintained - internally will be used. This is only useful in a limited number of - special cases, e.g. to describe the distance between remote points - of a single macromolecule, which cannot be split across periodic cell - boundaries, and for which the minimum-image distance might give the - wrong result because of a relatively small periodic cell.} - \item % \labelkey{colvar|distance|oneSiteTotalForce} \keydef @@ -1600,9 +1613,6 @@ (\bm{r}-\bm{r}_1)$. The vector should be written as three components separated by commas and enclosed in parentheses.} -\item % - \dupkey{forceNoPBC}{\texttt{distanceZ}}{colvar|distance|forceNoPBC}{\texttt{distance} component} - \item % \dupkey{oneSiteTotalForce}{\texttt{distanceZ}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1634,8 +1644,6 @@ \dupkey{ref2}{\texttt{distanceXY}}{colvar|distanceZ|ref2}{\texttt{distanceZ} component} \item % \dupkey{axis}{\texttt{distanceXY}}{colvar|distanceZ|axis}{\texttt{distanceZ} component} -\item % - \dupkey{forceNoPBC}{\texttt{distanceXY}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{distanceZ}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1645,19 +1653,17 @@ \cvsubsubsec{\texttt{distanceVec}: distance vector between two groups.}{sec:cvc_distanceVec} \labelkey{colvar|distanceVec} -The \texttt{distanceVec~\{...\}} block defines -a distance vector component, which accepts the same keywords as -the component \texttt{distance}: \texttt{group1}, \texttt{group2}, and -\texttt{forceNoPBC}. Its value is the 3-vector joining the centers -of mass of \texttt{group1} and \texttt{group2}. +The \texttt{distanceVec} component computes the 3-dimensional vector joining the centers of mass of +\texttt{group1} and \texttt{group2}. Its values are therefore multi-dimensional and are subject to +the restrictions listed in \ref{sec:cvc_non_scalar}. Moreover, when computing differences between +two different values of a \texttt{distanceVec} variable the minimum-image convention is assumed +(unless \refkey{forceNoPBC}{colvar|cvc|forceNoPBC} is enabled). \begin{cvcoptions} \item % \dupkey{group1}{\texttt{distanceVec}}{colvar|distance|group1}{\texttt{distance} component} \item % \simkey{group2}{\texttt{distanceVec}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{distanceVec}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{distanceVec}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1679,8 +1685,6 @@ \dupkey{group1}{\texttt{distanceDir}}{colvar|distance|group1}{\texttt{distance} component} \item % \simkey{group2}{\texttt{distanceDir}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{distanceDir}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{distanceDir}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1733,8 +1737,6 @@ \simkey{group2}{\texttt{angle}}{group1} \item % \simkey{group3}{\texttt{angle}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{angle}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{angle}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1756,8 +1758,6 @@ \simkey{group2}{\texttt{dipoleAngle}}{group1} \item % \simkey{group3}{\texttt{dipoleAngle}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{dipoleAngle}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{dipoleAngle}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -1783,8 +1783,6 @@ \simkey{group3}{\texttt{dihedral}}{group1} \item % \simkey{group4}{\texttt{dihedral}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{dihedral}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \item % \dupkey{oneSiteTotalForce}{\texttt{dihedral}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component} \end{cvcoptions} @@ -3039,8 +3037,6 @@ \dupkey{group1}{\texttt{distancePairs}}{colvar|distance|group1}{\texttt{distance} component} \item % \simkey{group2}{\texttt{distancePairs}}{group1} -\item % - \dupkey{forceNoPBC}{\texttt{distancePairs}}{colvar|distance|forceNoPBC}{\texttt{distance} component} \end{cvcoptions} This component returns a $N_{\mathrm{1}}\times{}N_{\mathrm{2}}$-dimensional vector of numbers, each ranging from $0$ to the largest possible distance within the chosen boundary conditions. @@ -3143,32 +3139,32 @@ The usage of \texttt{gzpath} and \texttt{gspath} is illustrated below: \begin{cvexampleinput} -colvar \{\\ -\-~~\# Progress along the path\\ -\-~~name gs\\ -\-~~\# The path is defined by 5 reference frames (from string-00.pdb to string-04.pdb)\\ -\-~~\# Use atomic coordinate from atoms 1, 2 and 3 to compute the path\\ -\-~~gspath \{\\ -\-~~~~atoms \{atomnumbers \{ 1 2 3 \}\}\\ -\-~~~~refPositionsFile1 string-00.pdb\\ -\-~~~~refPositionsFile2 string-01.pdb\\ -\-~~~~refPositionsFile3 string-02.pdb\\ -\-~~~~refPositionsFile4 string-03.pdb\\ -\-~~~~refPositionsFile5 string-04.pdb\\ +colvar~\{\\ +\-~~\#~Progress~along~the~path\\ +\-~~name~gs\\ +\-~~\#~The~path~is~defined~by~5~reference~frames~(from~string-00.pdb~to~string-04.pdb)\\ +\-~~\#~Use~atomic~coordinate~from~atoms~1,~2~and~3~to~compute~the~path\\ +\-~~gspath~\{\\ +\-~~~~atoms~\{atomnumbers~\{~1~2~3~\}\}\\ +\-~~~~refPositionsFile1~string-00.pdb\\ +\-~~~~refPositionsFile2~string-01.pdb\\ +\-~~~~refPositionsFile3~string-02.pdb\\ +\-~~~~refPositionsFile4~string-03.pdb\\ +\-~~~~refPositionsFile5~string-04.pdb\\ \-~~\}\\ \}\\ -\noindent\ttfamily colvar \{\\ -\-~~\# Distance from the path\\ -\-~~name gz\\ -\-~~\# The path is defined by 5 reference frames (from string-00.pdb to string-04.pdb)\\ -\-~~\# Use atomic coordinate from atoms 1, 2 and 3 to compute the path\\ -\-~~gzpath \{\\ -\-~~~~atoms \{atomnumbers \{ 1 2 3 \}\}\\ -\-~~~~refPositionsFile1 string-00.pdb\\ -\-~~~~refPositionsFile2 string-01.pdb\\ -\-~~~~refPositionsFile3 string-02.pdb\\ -\-~~~~refPositionsFile4 string-03.pdb\\ -\-~~~~refPositionsFile5 string-04.pdb\\ +\noindent\ttfamily~colvar~\{\\ +\-~~\#~Distance~from~the~path\\ +\-~~name~gz\\ +\-~~\#~The~path~is~defined~by~5~reference~frames~(from~string-00.pdb~to~string-04.pdb)\\ +\-~~\#~Use~atomic~coordinate~from~atoms~1,~2~and~3~to~compute~the~path\\ +\-~~gzpath~\{\\ +\-~~~~atoms~\{atomnumbers~\{~1~2~3~\}\}\\ +\-~~~~refPositionsFile1~string-00.pdb\\ +\-~~~~refPositionsFile2~string-01.pdb\\ +\-~~~~refPositionsFile3~string-02.pdb\\ +\-~~~~refPositionsFile4~string-03.pdb\\ +\-~~~~refPositionsFile5~string-04.pdb\\ \-~~\}\\ \} \end{cvexampleinput} @@ -3245,47 +3241,47 @@ The usage of \texttt{gzpathCV} and \texttt{gspathCV} is illustrated below: \begin{cvexampleinput} -\-colvar \{\\ -\-~~\# Progress along the path\\ -\-~~name gs\\ -\-~~\# Path defined by the CV space of two dihedral angles\\ -\-~~gspathCV \{\\ -\-~~~~pathFile ./path.txt\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 001\\ -\-~~~~~~group1 \{atomNumbers \{5\}\}\\ -\-~~~~~~group2 \{atomNumbers \{7\}\}\\ -\-~~~~~~group3 \{atomNumbers \{9\}\}\\ -\-~~~~~~group4 \{atomNumbers \{15\}\}\\ +\-colvar~\{\\ +\-~~\#~Progress~along~the~path\\ +\-~~name~gs\\ +\-~~\#~Path~defined~by~the~CV~space~of~two~dihedral~angles\\ +\-~~gspathCV~\{\\ +\-~~~~pathFile~./path.txt\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~001\\ +\-~~~~~~group1~\{atomNumbers~\{5\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{15\}\}\\ \-~~~~\}\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 002\\ -\-~~~~~~group1 \{atomNumbers \{7\}\}\\ -\-~~~~~~group2 \{atomNumbers \{9\}\}\\ -\-~~~~~~group3 \{atomNumbers \{15\}\}\\ -\-~~~~~~group4 \{atomNumbers \{17\}\}\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~002\\ +\-~~~~~~group1~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{15\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{17\}\}\\ \-~~~~\}\\ \-~~\}\\ \}\\ \\ -\-colvar \{\\ -\-~~\# Distance from the path\\ -\-~~name gz\\ -\-~~gzpathCV \{\\ -\-~~~~pathFile ./path.txt\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 001\\ -\-~~~~~~group1 \{atomNumbers \{5\}\}\\ -\-~~~~~~group2 \{atomNumbers \{7\}\}\\ -\-~~~~~~group3 \{atomNumbers \{9\}\}\\ -\-~~~~~~group4 \{atomNumbers \{15\}\}\\ +\-colvar~\{\\ +\-~~\#~Distance~from~the~path\\ +\-~~name~gz\\ +\-~~gzpathCV~\{\\ +\-~~~~pathFile~./path.txt\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~001\\ +\-~~~~~~group1~\{atomNumbers~\{5\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{15\}\}\\ \-~~~~\}\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 002\\ -\-~~~~~~group1 \{atomNumbers \{7\}\}\\ -\-~~~~~~group2 \{atomNumbers \{9\}\}\\ -\-~~~~~~group3 \{atomNumbers \{15\}\}\\ -\-~~~~~~group4 \{atomNumbers \{17\}\}\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~002\\ +\-~~~~~~group1~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{15\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{17\}\}\\ \-~~~~\}\\ \-~~\}\\ \} @@ -3366,51 +3362,51 @@ The usage of \texttt{azpathCV} and \texttt{aspathCV} is illustrated below: \begin{cvexampleinput} -colvar \{\\ -\-~~\# Progress along the path\\ -\-~~name as\\ -\-~~\# Path defined by the CV space of two dihedral angles\\ -\-~~aspathCV \{\\ -\-~~~~pathFile ./path.txt\\ -\-~~~~weights \{1.0 1.0\}\\ -\-~~~~lambda 0.005\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 001\\ -\-~~~~~~group1 \{atomNumbers \{5\}\}\\ -\-~~~~~~group2 \{atomNumbers \{7\}\}\\ -\-~~~~~~group3 \{atomNumbers \{9\}\}\\ -\-~~~~~~group4 \{atomNumbers \{15\}\}\\ +colvar~\{\\ +\-~~\#~Progress~along~the~path\\ +\-~~name~as\\ +\-~~\#~Path~defined~by~the~CV~space~of~two~dihedral~angles\\ +\-~~aspathCV~\{\\ +\-~~~~pathFile~./path.txt\\ +\-~~~~weights~\{1.0~1.0\}\\ +\-~~~~lambda~0.005\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~001\\ +\-~~~~~~group1~\{atomNumbers~\{5\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{15\}\}\\ \-~~~~\}\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 002\\ -\-~~~~~~group1 \{atomNumbers \{7\}\}\\ -\-~~~~~~group2 \{atomNumbers \{9\}\}\\ -\-~~~~~~group3 \{atomNumbers \{15\}\}\\ -\-~~~~~~group4 \{atomNumbers \{17\}\}\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~002\\ +\-~~~~~~group1~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{15\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{17\}\}\\ \-~~~~\}\\ \-~~\}\\ \}\\ \\ -colvar \{\\ -\-~~\# Distance from the path\\ -\-~~name az\\ -\-~~azpathCV \{\\ -\-~~~~pathFile ./path.txt\\ -\-~~~~weights \{1.0 1.0\}\\ -\-~~~~lambda 0.005\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 001\\ -\-~~~~~~group1 \{atomNumbers \{5\}\}\\ -\-~~~~~~group2 \{atomNumbers \{7\}\}\\ -\-~~~~~~group3 \{atomNumbers \{9\}\}\\ -\-~~~~~~group4 \{atomNumbers \{15\}\}\\ +colvar~\{\\ +\-~~\#~Distance~from~the~path\\ +\-~~name~az\\ +\-~~azpathCV~\{\\ +\-~~~~pathFile~./path.txt\\ +\-~~~~weights~\{1.0~1.0\}\\ +\-~~~~lambda~0.005\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~001\\ +\-~~~~~~group1~\{atomNumbers~\{5\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{15\}\}\\ \-~~~~\}\\ -\-~~~~dihedral \{\\ -\-~~~~~~name 002\\ -\-~~~~~~group1 \{atomNumbers \{7\}\}\\ -\-~~~~~~group2 \{atomNumbers \{9\}\}\\ -\-~~~~~~group3 \{atomNumbers \{15\}\}\\ -\-~~~~~~group4 \{atomNumbers \{17\}\}\\ +\-~~~~dihedral~\{\\ +\-~~~~~~name~002\\ +\-~~~~~~group1~\{atomNumbers~\{7\}\}\\ +\-~~~~~~group2~\{atomNumbers~\{9\}\}\\ +\-~~~~~~group3~\{atomNumbers~\{15\}\}\\ +\-~~~~~~group4~\{atomNumbers~\{17\}\}\\ \-~~~~\}\\ \-~~\}\\ \} @@ -3507,32 +3503,33 @@ The usage of \texttt{azpath} and \texttt{aspath} is illustrated below: \begin{cvexampleinput} -colvar \{\\ -\-~~\# Progress along the path\\ -\-~~name as\\ -\-~~\# The path is defined by 5 reference frames (from string-00.pdb to string-04.pdb)\\ -\-~~\# Use atomic coordinate from atoms 1, 2 and 3 to compute the path\\ -\-~~aspath \{\\ -\-~~~~atoms \{atomnumbers \{ 1 2 3 \}\}\\ -\-~~~~refPositionsFile1 string-00.pdb\\ -\-~~~~refPositionsFile2 string-01.pdb\\ -\-~~~~refPositionsFile3 string-02.pdb\\ -\-~~~~refPositionsFile4 string-03.pdb\\ -\-~~~~refPositionsFile5 string-04.pdb\\ +\-colvar~\{\\ +\-~~\#~Progress~along~the~path\\ +\-~~name~as\\ +\-~~\#~The~path~is~defined~by~5~reference~frames~(from~string-00.pdb~to~string-04.pdb)\\ +\-~~\#~Use~atomic~coordinate~from~atoms~1,~2~and~3~to~compute~the~path\\ +\-~~aspath~\{\\ +\-~~~~atoms~\{atomnumbers~\{~1~2~3~\}\}\\ +\-~~~~refPositionsFile1~string-00.pdb\\ +\-~~~~refPositionsFile2~string-01.pdb\\ +\-~~~~refPositionsFile3~string-02.pdb\\ +\-~~~~refPositionsFile4~string-03.pdb\\ +\-~~~~refPositionsFile5~string-04.pdb\\ \-~~\}\\ -\}\\ -\noindent\ttfamily colvar \{\\ -\-~~\# Distance from the path\\ -\-~~name az\\ -\-~~\# The path is defined by 5 reference frames (from string-00.pdb to string-04.pdb)\\ -\-~~\# Use atomic coordinate from atoms 1, 2 and 3 to compute the path\\ -\-~~azpath \{\\ -\-~~~~atoms \{atomnumbers \{ 1 2 3 \}\}\\ -\-~~~~refPositionsFile1 string-00.pdb\\ -\-~~~~refPositionsFile2 string-01.pdb\\ -\-~~~~refPositionsFile3 string-02.pdb\\ -\-~~~~refPositionsFile4 string-03.pdb\\ -\-~~~~refPositionsFile5 string-04.pdb\\ +\-\}\\ +\-\\ +\-colvar~\{\\ +\-~~\#~Distance~from~the~path\\ +\-~~name~az\\ +\-~~\#~The~path~is~defined~by~5~reference~frames~(from~string-00.pdb~to~string-04.pdb)\\ +\-~~\#~Use~atomic~coordinate~from~atoms~1,~2~and~3~to~compute~the~path\\ +\-~~azpath~\{\\ +\-~~~~atoms~\{atomnumbers~\{~1~2~3~\}\}\\ +\-~~~~refPositionsFile1~string-00.pdb\\ +\-~~~~refPositionsFile2~string-01.pdb\\ +\-~~~~refPositionsFile3~string-02.pdb\\ +\-~~~~refPositionsFile4~string-03.pdb\\ +\-~~~~refPositionsFile5~string-04.pdb\\ \-~~\}\\ \} \end{cvexampleinput} @@ -3873,42 +3870,40 @@ \cvsubsec{Periodic components}{sec:cvc_periodic} -The following components returns -real numbers that lie in a periodic interval: -\begin{itemize} -\item \texttt{dihedral}: torsional angle between four groups; -\item \texttt{spinAngle}: angle of rotation around a predefined axis - in the best-fit from a set of reference coordinates. -\end{itemize} -In certain conditions, \texttt{distanceZ} can also be periodic, namely -when periodic boundary conditions (PBCs) are defined in the simulation -and \texttt{distanceZ}'s axis is parallel to a unit cell vector. -In addition, a custom \cvscriptonly{or scripted} scalar colvar may be periodic -depending on its user-defined expression. It will only be treated as such by -the Colvars module if the period is specified using the \texttt{period} keyword, -while \texttt{wrapAround} is optional. +Certain components, such as \refkey{dihedral}{colvar|dihedral} or +\refkey{dihedral}{colvar|spinAngle}, compute angles that lie in a periodic interval between +$-180^\circ$ and $180^\circ$. When computing pairwise distances between values of those angles +(e.g.\ for the sake of computing restraint potentials, or sampling PMFs), periodicity is taken into +account by following the minimum-image convention. + +Additionally, several other components, such as \refkey{distanceZ}{colvar|distanceZ}, support +optional periodicity if this is provided in the configuration. -The following keywords can be used within periodic components\cvleptononly{, or within custom variables (\ref{sec:colvar_custom_function})}\cvscriptonly{, or wthin scripted variables \ref{sec:colvar_scripted}}). +The following keywords can be used within periodic components\cvleptononly{, or within custom + variables (\ref{sec:colvar_custom_function})}\cvscriptonly{, or wthin scripted variables + \ref{sec:colvar_scripted}}). \begin{itemize} \item % + \labelkey{colvar|cvc|period} \keydef {period}{% - \texttt{distanceZ}, custom colvars}{% + components with optional periodicity}{% Period of the component}{% positive decimal}{% - 0.0}{% - Setting this number enables the treatment of \texttt{distanceZ} as + 0.0, unless otherwise stated}{% + Setting this number enables the treatment of as a periodic component: by default, \texttt{distanceZ} is not considered periodic. The keyword is supported, but irrelevant within \texttt{dihedral} or \texttt{spinAngle}, because their period is always 360~degrees.} \item % + \labelkey{colvar|cvc|wrapAround} \keydef {wrapAround}{% - \texttt{distanceZ}, \texttt{dihedral}, \texttt{spinAngle}, custom colvars}{% + periodic components}{% Center of the wrapping interval for periodic variables}{% decimal}{% 0.0}{% @@ -3917,11 +3912,7 @@ This can be useful for convenience of output, or to set the walls for a \texttt{harmonicWalls} in an order that would not otherwise be allowed.} \end{itemize} -Internally, all differences between two values of a periodic colvar -follow the minimum image convention: they are calculated based on -the two periodic images that are closest to each other. - -\emph{Note: linear or polynomial combinations of periodic components (see \ref{sec:cvc_superp}) may become meaningless when components cross the periodic boundary. Use such combinations carefully: estimate the range of possible values of each component in a given simulation, and make use of \texttt{wrapAround} to limit this problem whenever possible.} +\emph{Note: using linear/polynomial combinations of periodic components (see \ref{sec:cvc_superp}), or other custom or scripted function may invalidate the periodicity. Use such combinations carefully: estimate the range of possible values of each component in a given simulation, and make use of \texttt{wrapAround} to limit this problem whenever possible.} \cvsubsec{Non-scalar components}{sec:cvc_non_scalar} @@ -4065,22 +4056,22 @@ A vector variable may be defined by specifying the \texttt{customFunction} parameter several times: each expression defines one scalar element of the vector colvar, as in this example: \begin{cvexampleinput} -\-colvar \{\\ -\-~~name custom\\ +\-colvar~\{\\ +\-~~name~custom\\ \\ -\-~~\# A 2-dimensional vector function of a scalar x and a 3-vector r\\ -\-~~customFunction cos(x) * (r1 + r2 + r3)\\ -\-~~customFunction sqrt(r1 * r2)\\ +\-~~\#~A~2-dimensional~vector~function~of~a~scalar~x~and~a~3-vector~r\\ +\-~~customFunction~cos(x)~*~(r1~+~r2~+~r3)\\ +\-~~customFunction~sqrt(r1~*~r2)\\ \\ -\-~~distance \{\\ -\-~~~~name x\\ -\-~~~~group1 \{ atomNumbers 1 \}\\ -\-~~~~group2 \{ atomNumbers 50 \}\\ +\-~~distance~\{\\ +\-~~~~name~x\\ +\-~~~~group1~\{~atomNumbers~1~\}\\ +\-~~~~group2~\{~atomNumbers~50~\}\\ \-~~\}\\ -\-~~distanceVec \{\\ -\-~~~~name r\\ -\-~~~~group1 \{ atomNumbers 10 11 12 \}\\ -\-~~~~group2 \{ atomNumbers 20 21 22 \}\\ +\-~~distanceVec~\{\\ +\-~~~~name~r\\ +\-~~~~group1~\{~atomNumbers~10~11~12~\}\\ +\-~~~~group2~\{~atomNumbers~~20~21~22~\}\\ \-~~\}\\ \} \end{cvexampleinput} @@ -4733,28 +4724,28 @@ \textbf{Note: }\emph{this is an unusually varied combination of selection keywords, demonstrating how they can be combined together: most simulations only use one of them.}\\ \begin{cvexampleinput} -\-atoms \{\\ +\-atoms~\{\\ \\ -\-~~\# add atoms 1 and 3 to this group (note: first atom in the system is 1)\\ -\-~~atomNumbers \{ \\ -\-~~~~1 3\\ +\-~~\#~add~atoms~1~and~3~to~this~group~(note:~first~atom~in~the~system~is~1)\\ +\-~~atomNumbers~\{~\\ +\-~~~~1~3\\ \-~~\}\\ \\ -\-~~\# add atoms starting from 20 up to and including 50\\ -\-~~atomNumbersRange 20-50\\ +\-~~\#~add~atoms~starting~from~20~up~to~and~including~50\\ +\-~~atomNumbersRange~~20-50\\ \cvnamebasedonly{% -\-~~\# add all the atoms with occupancy 2 in the file atoms.pdb\\ -\-~~atomsFile atoms.pdb\\ -\-~~atomsCol O\\ -\-~~atomsColValue 2.0\\ +\-~~\#~add~all~the~atoms~with~occupancy~2~in~the~file~atoms.pdb\\ +\-~~atomsFile~~~~~~~~~~~~~atoms.pdb\\ +\-~~atomsCol~~~~~~~~~~~~~~O\\ +\-~~atomsColValue~~~~~~~~~2.0\\ \\ -\-~~\# add all the C-alphas within residues 11 to 20 of segments "PR1" and "PR2"\\ -\-~~psfSegID PR1 PR2\\ -\-~~atomNameResidueRange CA 11-20\\ -\-~~atomNameResidueRange CA 11-20\\ +\-~~\#~add~all~the~C-alphas~within~residues~11~to~20~of~segments~"PR1"~and~"PR2"\\ +\-~~psfSegID~~~~~~~~~~~~~~PR1~PR2\\ +\-~~atomNameResidueRange~~CA~11-20\\ +\-~~atomNameResidueRange~~CA~11-20\\ } -\-~~\# add index group (requires a .ndx file to be provided globally)\\ -\-~~indexGroup Water\\ +\-~~\#~add~index~group~(requires~a~.ndx~file~to~be~provided~globally)\\ +\-~~indexGroup~Water\\ \} \end{cvexampleinput} @@ -5002,22 +4993,22 @@ To define a DBC coordinate, this atom group would be used within an \refkey{rmsd}{sec:cvc_rmsd} function. \begin{cvexampleinput} -\-\# Example: defining a group "atoms" (the ligand) whose coordinates are expressed \\ -\# in a roto-translated frame of reference defined by a second group (the receptor)\\ -atoms \{\\ +\-\#~Example:~defining~a~group~"atoms"~(the~ligand)~whose~coordinates~are~expressed~\\ +\#~in~a~roto-translated~frame~of~reference~defined~by~a~second~group~(the~receptor)\\ +atoms~\{\\ \\ -\-~~atomNumbers 1 2 3 4 5 6 7 \# atoms of the ligand (1-based)\\ +\-~~atomNumbers~1~2~3~4~5~6~7~\#~atoms~of~the~ligand~(1-based)\\ \\ -\-~~centerToReference yes\\ -\-~~rotateToReference yes\\ -\-~~fittingGroup \{\\ -\-~~~~\# define the frame by fitting alpha carbon atoms\\ -\-~~~~\# in 2 protein segments close to the site\\ -\-~~~~psfSegID PROT PROT\\ -\-~~~~atomNameResidueRange CA 1-40\\ -\-~~~~atomNameResidueRange CA 59-100\\ -\-~~\} \\ -\-~~refPositionsFile all.pdb \# can be the entire system\\ +\-~~centerToReference~yes\\ +\-~~rotateToReference~yes\\ +\-~~fittingGroup~\{\\ +\-~~~~\#~define~the~frame~by~fitting~alpha~carbon~atoms\\ +\-~~~~\#~in~2~protein~segments~close~to~the~site\\ +\-~~~~psfSegID~~~~~~~~~~~~~~PROT~PROT\\ +\-~~~~atomNameResidueRange~~CA~1-40\\ +\-~~~~atomNameResidueRange~~CA~59-100\\ +\-~~\}~\\ +\-~~refPositionsFile~all.pdb~~\#~can~be~the~entire~system\\ \} \end{cvexampleinput} } @@ -5133,7 +5124,6 @@ A biasing or analysis method can be applied to existing collective variables by using the following configuration: \begin{cvexampleinput} - \-$<$biastype$>$~\{\\ \-~~name~$<$name$>$ \\ \-~~colvars~$<$xi1$>$~$<$xi2$>$~...\\ @@ -5970,25 +5960,25 @@ \noindent\textbf{Example:} Using harmonic walls to protect the grid's boundaries. \begin{cvexampleinput} -\-colvar \{\\ -\-\-~~name r\\ -\-\-~~distance \{ ... \}\\ -\-\-~~upperBoundary 15.0\\ -\-\-~~width 0.2\\ +\-colvar~\{\\ +\-\-~~name~r\\ +\-\-~~distance~\{~...~\}\\ +\-\-~~upperBoundary~15.0\\ +\-\-~~width~0.2\\ \-\}\\ \\ -\-metadynamics \{\\ -\-\-~~name meta\_r\\ -\-\-~~colvars r\\ -\-\-~~hillWeight 0.001\\ -\-\-~~hillWidth 2.0\\ +\-metadynamics~\{\\ +\-\-~~name~meta\_r\\ +\-\-~~colvars~r\\ +\-\-~~hillWeight~0.001\\ +\-\-~~hillWidth~2.0\\ \-\}\\ \\ -\-harmonicWalls \{\\ -\-\-~~name wall\_r\\ -\-\-~~colvars r\\ -\-\-~~upperWalls 13.0\\ -\-\-~~upperWallConstant 2.0\\ +\-harmonicWalls~\{\\ +\-\-~~name~wall\_r\\ +\-\-~~colvars~r\\ +\-\-~~upperWalls~13.0\\ +\-\-~~upperWallConstant~2.0\\ \-\} \end{cvexampleinput} @@ -6271,24 +6261,24 @@ \noindent\textbf{Example:} EBmetaD configuration for a single variable. \begin{cvexampleinput} -\-colvar \{\\ -\-\-~~name r \\ -\-\-~~distance \{\\ -\-\-~~~~group1 \{ atomNumbers 991 992 \}\\ -\-\-~~~~group2 \{ atomNumbers 1762 1763 \}\\ +\-colvar~\{\\ +\-\-~~name~r~\\ +\-\-~~distance~\{\\ +\-\-~~~~group1~\{~atomNumbers~991~992~\}\\ +\-\-~~~~group2~\{~atomNumbers~1762~1763~\}\\ \-\-~~\}\\ -\-\-~~upperBoundary 100.0 \\ -\-\-~~width 0.1 \\ +\-\-~~upperBoundary~~100.0~\\ +\-\-~~width~~~~~~~~~~~~0.1~\\ \-\}\\ \\ -\-metadynamics \{\\ -\-\-~~name ebmeta\\ -\-\-~~colvars r\\ -\-\-~~hillWeight 0.01\\ -\-\-~~hillWidth 3.0\\ -\-\-~~ebMeta on\\ -\-\-~~targetDistFile targetdist1.dat\\ -\-\-~~ebMetaEquilSteps 500000\\ +\-metadynamics~\{\\ +\-\-~~name~~~~~~~~~~~~~~ebmeta\\ +\-\-~~colvars~~~~~~~~~~~r\\ +\-\-~~hillWeight~~~~~~~~0.01\\ +\-\-~~hillWidth~~~~~~~~~3.0\\ +\-\-~~ebMeta~~~~~~~~~~~~on\\ +\-\-~~targetDistFile~~~~targetdist1.dat\\ +\-\-~~ebMetaEquilSteps~~500000\\ \-\} \end{cvexampleinput} @@ -6363,16 +6353,16 @@ \noindent\textbf{Example:} Multiple-walker metadynamics with file-based communication.\\ \begin{cvexampleinput} -\-metadynamics \{\\ -\-\-~~name mymtd\\ -\-\-~~colvars x\\ -\-\-~~hillWeight 0.001\\ -\-\-~~newHillFrequency 1000\\ -\-\-~~hillWidth 3.0\\ +\-metadynamics~\{\\ +\-\-~~name~mymtd\\ +\-\-~~colvars~x\\ +\-\-~~hillWeight~0.001\\ +\-\-~~newHillFrequency~1000\\ +\-\-~~hillWidth~3.0\\ \-\-~~\\ -\-\-~~multipleReplicas on\\ -\-\-~~replicasRegistry /shared-folder/mymtd-replicas.txt\\ -\-\-~~replicaUpdateFrequency 50000 \# Best if larger than newHillFrequency\\ +\-\-~~multipleReplicas~~~~~~~on\\ +\-\-~~replicasRegistry~~~~~~~/shared-folder/mymtd-replicas.txt\\ +\-\-~~replicaUpdateFrequency~50000~~\#~Best~if~larger~than~newHillFrequency\\ \-\} \end{cvexampleinput} @@ -6827,24 +6817,24 @@ \noindent\textbf{Example 1:} harmonic walls for one variable with two different force constants. \begin{cvexampleinput} -\-harmonicWalls \{\\ -\-\-~~name mywalls\\ -\-\-~~colvars dist\\ -\-\-~~lowerWalls 22.0 \\ -\-\-~~upperWalls 38.0 \\ -\-\-~~lowerWallConstant 2.0 \\ -\-\-~~upperWallConstant 10.0 \\ +\-harmonicWalls~\{\\ +\-\-~~name~~mywalls\\ +\-\-~~colvars~~~~dist\\ +\-\-~~lowerWalls~~22.0~\\ +\-\-~~upperWalls~~38.0~\\ +\-\-~~lowerWallConstant~~2.0~\\ +\-\-~~upperWallConstant~10.0~\\ \-\} \end{cvexampleinput} \noindent\textbf{Example 2:} harmonic walls for two variables with a single force constant. \begin{cvexampleinput} -\-harmonicWalls \{\\ -\-\-~~name mywalls\\ -\-\-~~colvars phi psi\\ -\-\-~~lowerWalls -180.0 0.0\\ -\-\-~~upperWalls 0.0 180.0\\ -\-\-~~forceConstant 5.0 \\ +\-harmonicWalls~\{\\ +\-\-~~name~~mywalls\\ +\-\-~~colvars~~~~~~~phi~~~~psi\\ +\-\-~~lowerWalls~-180.0~~~~0.0\\ +\-\-~~upperWalls~~~~0.0~~180.0\\ +\-\-~~forceConstant~5.0~\\ \-\} \end{cvexampleinput} @@ -7482,3 +7472,7 @@ \url{https://colvars.github.io/README-c++11.html} \end{itemize} +% Emacs +% Local Variables: +% fill-column: 100 +% End: diff --git a/doc/colvars-refman.tex b/doc/colvars-refman.tex index 001a46d00..245a33b0e 100644 --- a/doc/colvars-refman.tex +++ b/doc/colvars-refman.tex @@ -232,13 +232,14 @@ \newcommand{\bm}[1]{{\mbox{\boldmath{$#1$}}}} \newenvironment{cvcoptions}{% - %% TODO Uncomment when all itemizations environment are redefined - %% \begin{mdframed}[linecolor=cvexample-border-color,linewidth=2pt] - \noindent\textbf{List of keywords} (see also \ref{sec:cvc_superp} for additional options): - \begin{itemize}}{% + \noindent\textbf{List of keywords} (see also \ref{sec:cvc_pbcs}, \ref{sec:cvc_common}, \ref{sec:cvc_periodic} and \ref{sec:cvc_superp} for additional options): + \begin{itemize} + }{% + \item % + \colorbox{cvkeyword-background-color}{\texttt{forceNoPBC}} --- \emph{Use absolute rather than minimum-image distances?} (see \hyperlink{colvar|cvc|forceNoPBC}{description})\\ + Default: \texttt{no} \end{itemize} - %% \end{mdframed} - } +} \newcommand{\pagetoplink}{\ifdefined\HCode \begin{flushright} diff --git a/doc/extract_code_refs.py b/doc/extract_code_refs.py index 3f3a7f726..5de841ff3 100755 --- a/doc/extract_code_refs.py +++ b/doc/extract_code_refs.py @@ -8,7 +8,7 @@ def gen_cplusplus_entry(entry, key, url): - key_year_sep = re.search('\d', key).start() + key_year_sep = re.search('\\d', key).start() key_author = key[0:key_year_sep] key_year = key[key_year_sep:] result = """ diff --git a/namd/tests/library/033_apathCV_distanceVec_derivatives/apath.in b/namd/tests/library/033_apathCV_distanceVec_derivatives/apath.in index 626251f74..820e62020 100644 --- a/namd/tests/library/033_apathCV_distanceVec_derivatives/apath.in +++ b/namd/tests/library/033_apathCV_distanceVec_derivatives/apath.in @@ -9,7 +9,6 @@ colvar { distanceVec { name 002 componentCoeff 0.6 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -19,7 +18,6 @@ colvar { } distanceVec { name 005 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -29,7 +27,6 @@ colvar { } distanceVec { name 006 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -39,7 +36,6 @@ colvar { } distanceVec { name 007 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -49,7 +45,6 @@ colvar { } distanceVec { name 009 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -59,7 +54,6 @@ colvar { } distanceVec { name 011 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -69,7 +63,6 @@ colvar { } distanceVec { name 015 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -79,7 +72,6 @@ colvar { } distanceVec { name 016 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -89,7 +81,6 @@ colvar { } distanceVec { name 017 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -99,7 +90,6 @@ colvar { } distanceVec { name 019 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -109,7 +99,6 @@ colvar { } distanceVec { name 021 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -119,7 +108,6 @@ colvar { } distanceVec { name 025 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -129,7 +117,6 @@ colvar { } distanceVec { name 026 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -139,7 +126,6 @@ colvar { } distanceVec { name 027 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -149,7 +135,6 @@ colvar { } distanceVec { name 029 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -159,7 +144,6 @@ colvar { } distanceVec { name 031 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -169,7 +153,6 @@ colvar { } distanceVec { name 035 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -179,7 +162,6 @@ colvar { } distanceVec { name 036 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -189,7 +171,6 @@ colvar { } distanceVec { name 037 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -199,7 +180,6 @@ colvar { } distanceVec { name 039 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -219,7 +199,6 @@ colvar { debugGradients on distanceVec { name 002 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -229,7 +208,6 @@ colvar { } distanceVec { name 005 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -239,7 +217,6 @@ colvar { } distanceVec { name 006 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -249,7 +226,6 @@ colvar { } distanceVec { name 007 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -259,7 +235,6 @@ colvar { } distanceVec { name 009 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -269,7 +244,6 @@ colvar { } distanceVec { name 011 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -279,7 +253,6 @@ colvar { } distanceVec { name 015 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -289,7 +262,6 @@ colvar { } distanceVec { name 016 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -299,7 +271,6 @@ colvar { } distanceVec { name 017 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -309,7 +280,6 @@ colvar { } distanceVec { name 019 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -319,7 +289,6 @@ colvar { } distanceVec { name 021 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -329,7 +298,6 @@ colvar { } distanceVec { name 025 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -339,7 +307,6 @@ colvar { } distanceVec { name 026 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -349,7 +316,6 @@ colvar { } distanceVec { name 027 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -359,7 +325,6 @@ colvar { } distanceVec { name 029 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -369,7 +334,6 @@ colvar { } distanceVec { name 031 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -379,7 +343,6 @@ colvar { } distanceVec { name 035 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -389,7 +352,6 @@ colvar { } distanceVec { name 036 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -399,7 +361,6 @@ colvar { } distanceVec { name 037 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } @@ -409,7 +370,6 @@ colvar { } distanceVec { name 039 - forceNoPBC yes group1 { dummyAtom (0.00, 0.00, 0.00) } diff --git a/src/colvar.cpp b/src/colvar.cpp index cc0d0e3f2..37911cdb7 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -22,10 +22,8 @@ #include "colvars_memstream.h" - -std::map> - colvar::global_cvc_map = - std::map>(); +std::map> colvar::global_cvc_map = + std::map>(); std::map colvar::global_cvc_desc_map = std::map(); @@ -139,7 +137,14 @@ int colvar::init(std::string const &conf) // Sort array of cvcs based on their names // Note: default CVC names are in input order for same type of CVC - std::sort(cvcs.begin(), cvcs.end(), colvar::compare_cvc); + std::sort(cvcs.begin(), cvcs.end(), + [](std::shared_ptr const &cvc1, + std::shared_ptr const &cvc2) -> bool { + if (cvc1 && cvc2) { + return cvc1->name < cvc2->name; + } + return false; + }); if(cvcs.size() > 1) { cvm::log("Sorted list of components for this scripted colvar:\n"); @@ -194,9 +199,9 @@ int colvar::init(std::string const &conf) if ((cvcs[i])->sup_np < 0) { cvm::log("Warning: you chose a negative exponent in the combination; " - "if you apply forces, the simulation may become unstable " - "when the component \""+ - (cvcs[i])->function_type+"\" approaches zero.\n"); + "if you apply forces, the simulation may become unstable " + "when the component \""+ + (cvcs[i])->function_type()+"\" approaches zero.\n"); } } } @@ -303,7 +308,7 @@ int colvar::init(std::string const &conf) error_code |= init_grid_parameters(conf); // Detect if we have a single component that is an alchemical lambda - if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alchLambda") { + if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type() == "alchLambda") { enable(f_cv_external); } @@ -755,8 +760,8 @@ template void colvar::add_component_type(char const *def_description, char const *def_config_key) { if (global_cvc_map.count(def_config_key) == 0) { - global_cvc_map[def_config_key] = [](const std::string &cvc_conf) { - return new def_class_name(cvc_conf); + global_cvc_map[def_config_key] = []() { + return new def_class_name(); }; global_cvc_desc_map[def_config_key] = std::string(def_description); } @@ -767,6 +772,7 @@ int colvar::init_components_type(const std::string& conf, const char* def_config size_t def_count = 0; std::string def_conf = ""; size_t pos = 0; + int error_code = COLVARS_OK; while ( this->key_lookup(conf, def_config_key, &def_conf, @@ -776,39 +782,29 @@ int colvar::init_components_type(const std::string& conf, const char* def_config "a new \""+std::string(def_config_key)+"\" component"+ (cvm::debug() ? ", with configuration:\n"+def_conf : ".\n")); - cvc *cvcp = global_cvc_map[def_config_key](def_conf); - cvm::increase_depth(); - if (cvcp) { - int error_code = cvcp->init_code; - cvcs.push_back(cvcp); - error_code |= cvcp->set_function_type(def_config_key); - if (error_code == COLVARS_OK) { - error_code |= cvcp->check_keywords(def_conf, def_config_key); - } - if (error_code != COLVARS_OK) { - cvm::decrease_depth(); - return cvm::error("Error: in setting up component \"" + std::string(def_config_key) + - "\".\n", - COLVARS_INPUT_ERROR); - } - } else { - cvm::decrease_depth(); - return cvm::error("Error: in allocating component \"" + std::string(def_config_key) + "\".\n", + cvc *cvcp = global_cvc_map[def_config_key](); + if (!cvcp) { + return cvm::error("Error: in creating object of type \"" + std::string(def_config_key) + + "\".\n", COLVARS_MEMORY_ERROR); } + cvcs.push_back(std::shared_ptr(cvcp)); - if ((cvcp->period != 0.0) || (cvcp->wrap_center != 0.0)) { - if (!cvcp->is_enabled(f_cvc_periodic)) { - cvm::decrease_depth(); - return cvm::error("Error: invalid use of period and/or " - "wrapAround in a \"" + - std::string(def_config_key) + "\" component.\n" + - "Period: " + cvm::to_str(cvcp->period) + - " wrapAround: " + cvm::to_str(cvcp->wrap_center), - COLVARS_INPUT_ERROR); - } + cvm::increase_depth(); + int error_code_this = cvcp->init(def_conf); + if (error_code_this == COLVARS_OK) { + // Checking for invalid keywords only if the parsing was successful, otherwise any + // early-returns due to errors would raise false positives + error_code_this |= cvcp->check_keywords(def_conf, def_config_key); + } + cvm::decrease_depth(); + if (error_code_this != COLVARS_OK) { + error_code |= + cvm::error("Error: in setting up component \"" + std::string(def_config_key) + "\".\n", + COLVARS_INPUT_ERROR); } + // Set default name if it doesn't have one if ( ! cvcs.back()->name.size()) { std::ostringstream s; s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count; @@ -822,15 +818,13 @@ int colvar::init_components_type(const std::string& conf, const char* def_config (cvm::debug() ? ", named \"" + cvcs.back()->name + "\"" : "") + ".\n"); } - cvm::decrease_depth(); - def_conf = ""; if (cvm::debug()) { cvm::log("Parsed " + cvm::to_str(cvcs.size()) + " components at this time.\n"); } } - return COLVARS_OK; + return error_code; } @@ -944,7 +938,7 @@ int colvar::init_components(std::string const &conf) if (error_code == COLVARS_OK) { // Store list of children cvcs for dependency checking purposes for (i = 0; i < cvcs.size(); i++) { - add_child(cvcs[i]); + add_child(cvcs[i].get()); } // By default all CVCs are active at the start n_active_cvcs = cvcs.size(); @@ -1284,14 +1278,10 @@ colvar::~colvar() // for dependency purposes remove_all_children(); - for (std::vector::reverse_iterator ci = cvcs.rbegin(); - ci != cvcs.rend(); - ++ci) { - // clear all children of this cvc (i.e. its atom groups) - // because the cvc base class destructor can't do it early enough - // and we don't want to have each cvc derived class do it separately + for (auto ci = cvcs.rbegin(); ci != cvcs.rend(); ++ci) { + // Clear all children of this cvc (i.e. its atom groups), because the cvc base class destructor + // can't do it early enough and we don't want to have each cvc derived class do it separately (*ci)->remove_all_children(); - delete *ci; } cvcs.clear(); diff --git a/src/colvar.h b/src/colvar.h index 0cc4208e1..443e1e4bd 100644 --- a/src/colvar.h +++ b/src/colvar.h @@ -14,6 +14,7 @@ #include #include #include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -617,7 +618,6 @@ class colvar : public colvarparse, public colvardeps { class dihedPC; class alch_lambda; class alch_Flambda; - class componentDisabled; class CartesianBasedPath; class aspath; class azpath; @@ -645,8 +645,7 @@ class colvar : public colvarparse, public colvardeps { class map_total; /// A global mapping of cvc names to the cvc constructors - static const std::map> & - get_global_cvc_map() + static const std::map> &get_global_cvc_map() { return global_cvc_map; } @@ -656,8 +655,8 @@ class colvar : public colvarparse, public colvardeps { protected: - /// \brief Array of \link colvar::cvc \endlink objects - std::vector cvcs; + /// Array of components objects + std::vector> cvcs; /// \brief Flags to enable or disable cvcs at next colvar evaluation std::vector cvc_flags; @@ -689,8 +688,7 @@ class colvar : public colvarparse, public colvardeps { #endif /// A global mapping of cvc names to the cvc constructors - static std::map> - global_cvc_map; + static std::map> global_cvc_map; /// A global mapping of cvc names to the corresponding descriptions static std::map global_cvc_desc_map; diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 4c7ab09dc..e6729f43a 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -19,63 +19,40 @@ colvar::cvc::cvc() { description = "uninitialized colvar component"; - b_try_scalable = true; - sup_coeff = 1.0; - sup_np = 1; - period = 0.0; - wrap_center = 0.0; - width = 0.0; cvc::init_dependencies(); } -colvar::cvc::cvc(std::string const &conf) -{ - description = "uninitialized colvar component"; - b_try_scalable = true; - sup_coeff = 1.0; - sup_np = 1; - period = 0.0; - wrap_center = 0.0; - width = 0.0; - init_dependencies(); - colvar::cvc::init(conf); -} - - int colvar::cvc::update_description() { if (name.size() > 0) { - description = "cvc " + name; + description = "cvc \"" + name + "\""; } else { description = "unnamed cvc"; } - if (function_type.size() > 0) { - description += " of type \"" + function_type + "\""; - } else { - description += " of unset type"; - } + description += " of type \"" + function_type() + "\""; return COLVARS_OK; } -int colvar::cvc::set_function_type(std::string const &type) +std::string colvar::cvc::function_type() const { - function_type = type; - if (function_types.size() == 0) { - function_types.push_back(function_type); - } else { - if (function_types.back() != function_type) { - function_types.push_back(function_type); - } + if (function_types.empty()) { + return "unset"; } - update_description(); + return function_types.back(); +} + +int colvar::cvc::set_function_type(std::string const &type) +{ + function_types.push_back(type); + update_description(); + cvm::main()->cite_feature(function_types[0]+" colvar component"); for (size_t i = function_types.size()-1; i > 0; i--) { cvm::main()->cite_feature(function_types[i]+" colvar component"+ " (derived from "+function_types[i-1]+")"); } - cvm::main()->cite_feature(function_types[0]+" colvar component"); return COLVARS_OK; } @@ -85,6 +62,8 @@ int colvar::cvc::init(std::string const &conf) if (cvm::debug()) cvm::log("Initializing cvc base object.\n"); + int error_code = COLVARS_OK; + std::string const old_name(name); if (name.size() > 0) { @@ -93,9 +72,9 @@ int colvar::cvc::init(std::string const &conf) if (get_keyval(conf, "name", name, name)) { if ((name != old_name) && (old_name.size() > 0)) { - cvm::error("Error: cannot rename component \""+old_name+ - "\" after initialization (new name = \""+name+"\")", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: cannot rename component \"" + old_name + + "\" after initialization (new name = \"" + name + "\")", + COLVARS_INPUT_ERROR); name = old_name; } } @@ -116,6 +95,24 @@ int colvar::cvc::init(std::string const &conf) register_param("period", reinterpret_cast(&period)); register_param("wrapAround", reinterpret_cast(&wrap_center)); + if (period != 0.0) { + if (!is_available(f_cvc_periodic)) { + error_code |= + cvm::error("Error: invalid use of period and/or " + "wrapAround in a \"" + + function_type() + "\" component.\n" + "Period: " + cvm::to_str(period) + + " wrapAround: " + cvm::to_str(wrap_center), + COLVARS_INPUT_ERROR); + } else { + enable(f_cvc_periodic); + } + } + + if ((wrap_center != 0.0) && !is_enabled(f_cvc_periodic)) { + error_code |= cvm::error("Error: wrapAround was defined for a non-periodic component.\n", + COLVARS_INPUT_ERROR); + } + get_keyval_feature(this, conf, "debugGradients", f_cvc_debug_gradient, false, parse_silent); @@ -133,7 +130,7 @@ int colvar::cvc::init(std::string const &conf) if (cvm::debug()) cvm::log("Done initializing cvc base object.\n"); - return cvm::get_error(); + return error_code; } @@ -171,12 +168,13 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, char const *group_key, bool optional) { - int &error_code = init_code; + int error_code = COLVARS_OK; cvm::atom_group *group = nullptr; std::string group_conf; if (key_lookup(conf, group_key, &group_conf)) { + group = new cvm::atom_group(group_key); if (b_try_scalable) { @@ -198,6 +196,8 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, COLVARS_INPUT_ERROR); delete group; group = nullptr; + // Silence unused variable warning; TODO stop returning a pointer + (void) error_code; return group; } @@ -224,6 +224,9 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, } } + // Silence unused variable warning; TODO stop returning a pointer + (void) error_code; + return group; } @@ -250,6 +253,8 @@ int colvar::cvc::init_dependencies() { init_feature(f_cvc_upper_boundary, "defined_upper_boundary", f_type_static); + init_feature(f_cvc_explicit_atom_groups, "explicit_atom_groups", f_type_static); + init_feature(f_cvc_gradient, "gradient", f_type_dynamic); init_feature(f_cvc_explicit_gradient, "explicit_gradient", f_type_static); @@ -286,6 +291,7 @@ int colvar::cvc::init_dependencies() { init_feature(f_cvc_collect_atom_ids, "collect_atom_ids", f_type_dynamic); require_feature_children(f_cvc_collect_atom_ids, f_ag_collect_atom_ids); + require_feature_self(f_cvc_collect_atom_ids, f_cvc_explicit_atom_groups); // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group // require_feature_children(f_cvc_scalable, f_ag_scalable); @@ -314,6 +320,8 @@ int colvar::cvc::init_dependencies() { feature_states[f_cvc_gradient].available = true; feature_states[f_cvc_collect_atom_ids].available = true; + feature_states[f_cvc_periodic].available = false; + // CVCs are enabled from the start - get disabled based on flags enable(f_cvc_active); @@ -371,6 +379,7 @@ void colvar::cvc::init_as_angle() void colvar::cvc::init_as_periodic_angle() { x.type(colvarvalue::type_scalar); + provide(f_cvc_periodic); enable(f_cvc_periodic); period = 360.0; init_scalar_boundaries(-180.0, 180.0); @@ -394,6 +403,7 @@ void colvar::cvc::register_atom_group(cvm::atom_group *ag) { atom_groups.push_back(ag); add_child(ag); + enable(f_cvc_explicit_atom_groups); } @@ -433,29 +443,21 @@ int colvar::cvc::set_param(std::string const ¶m_name, void colvar::cvc::read_data() { - size_t ig; - for (ig = 0; ig < atom_groups.size(); ig++) { - cvm::atom_group &atoms = *(atom_groups[ig]); - atoms.reset_atoms_data(); - atoms.read_positions(); - atoms.calc_required_properties(); - // each atom group will take care of its own fitting_group, if defined + if (is_enabled(f_cvc_explicit_atom_groups)) { + for (auto agi = atom_groups.begin(); agi != atom_groups.end(); agi++) { + cvm::atom_group &atoms = *(*agi); + atoms.reset_atoms_data(); + atoms.read_positions(); + atoms.calc_required_properties(); + // each atom group will take care of its own fitting_group, if defined + } } - -//// Don't try to get atom velocities, as no back-end currently implements it -// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) { -// for (i = 0; i < cvcs.size(); i++) { -// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { -// cvcs[i]->atom_groups[ig]->read_velocities(); -// } -// } -// } } -std::vector > colvar::cvc::get_atom_lists() +std::vector> colvar::cvc::get_atom_lists() { - std::vector > lists; + std::vector> lists; std::vector::iterator agi = atom_groups.begin(); for ( ; agi != atom_groups.end(); ++agi) { @@ -516,7 +518,7 @@ void colvar::cvc::collect_gradients(std::vector const &atom_ids, std::vecto void colvar::cvc::calc_force_invgrads() { cvm::error("Error: calculation of inverse gradients is not implemented " - "for colvar components of type \""+function_type+"\".\n", + "for colvar components of type \""+function_type()+"\".\n", COLVARS_NOT_IMPLEMENTED); } @@ -524,7 +526,7 @@ void colvar::cvc::calc_force_invgrads() void colvar::cvc::calc_Jacobian_derivative() { cvm::error("Error: calculation of inverse gradients is not implemented " - "for colvar components of type \""+function_type+"\".\n", + "for colvar components of type \""+function_type()+"\".\n", COLVARS_NOT_IMPLEMENTED); } @@ -537,6 +539,18 @@ void colvar::cvc::calc_fit_gradients() } +void colvar::cvc::apply_force(colvarvalue const &cvforce) +{ + if (is_enabled(f_cvc_explicit_atom_groups)) { + for (auto agi = atom_groups.begin(); agi != atom_groups.end(); agi++) { + if (!(*agi)->noforce) { + (*agi)->apply_colvar_force(cvforce); + } + } + } +} + + void colvar::cvc::debug_gradients() { // this function should work for any scalar cvc: @@ -656,34 +670,43 @@ void colvar::cvc::debug_gradients() } -cvm::real colvar::cvc::dist2(colvarvalue const &x1, - colvarvalue const &x2) const +cvm::real colvar::cvc::dist2(colvarvalue const &x1, colvarvalue const &x2) const { - return x1.dist2(x2); + cvm::real diff = x1.real_value - x2.real_value; + if (is_enabled(f_cvc_periodic)) { + cvm::real const shift = cvm::floor(diff / period + 0.5); + diff -= shift * period; + } + return diff * diff; } -colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return x1.dist2_grad(x2); + cvm::real diff = x1.real_value - x2.real_value; + if (is_enabled(f_cvc_periodic)) { + cvm::real const shift = cvm::floor(diff / period + 0.5); + diff -= shift * period; + } + return 2.0 * diff; } -colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return x2.dist2_grad(x1); + return cvc::dist2_lgrad(x1, x2); } -void colvar::cvc::wrap(colvarvalue & /* x_unwrapped */) const +void colvar::cvc::wrap(colvarvalue &x_unwrapped) const { - return; + if (is_enabled(f_cvc_periodic)) { + cvm::real const shift = cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5); + x_unwrapped.real_value -= shift * period; + } } - // Static members std::vector colvar::cvc::cvc_features; diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 1b0c631f4..334fdc1f6 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -74,40 +74,29 @@ class colvar::cvc /// cvc instance when debugging) std::string name; - /// \brief Description of the type of collective variable - /// - /// Normally this string is set by the parent \link colvar \endlink - /// object within its constructor, when all \link colvar::cvc \endlink - /// objects are initialized; therefore the main "config string" - /// constructor does not need to define it. If a \link colvar::cvc - /// \endlink is initialized and/or a different constructor is used, - /// this variable definition should be set within the constructor. - std::string function_type; + /// String identifier for the type of collective variable + std::string function_type() const; /// Keyword used in the input to denote this CVC std::string config_key; /// \brief Coefficient in the polynomial combination (default: 1.0) - cvm::real sup_coeff; + cvm::real sup_coeff = 1.0; + /// \brief Exponent in the polynomial combination (default: 1) - int sup_np; + int sup_np = 1; /// \brief Period of the values of this CVC (default: 0.0, non periodic) - cvm::real period; + cvm::real period = 0.0; /// \brief If the component is periodic, wrap around this value (default: 0.0) - cvm::real wrap_center; - - /// \brief Constructor - /// - /// Calls the init() function of the class - cvc(std::string const &conf); + cvm::real wrap_center = 0.0; - /// Current initialization state; TODO remove this when using init() after default constructor - int init_code = COLVARS_OK; + /// Constructor + cvc(); - /// Set the value of \link function_type \endlink and its dependencies - int set_function_type(std::string const &type); + /// Destructor + virtual ~cvc(); /// An init function should be defined for every class inheriting from cvc /// \param conf Contents of the configuration file pertaining to this \link @@ -117,38 +106,20 @@ class colvar::cvc /// \brief Initialize dependency tree virtual int init_dependencies(); - /// \brief Within the constructor, make a group parse its own - /// options from the provided configuration string - /// Returns reference to new group - cvm::atom_group *parse_group(std::string const &conf, - char const *group_key, - bool optional = false); - - /// \brief Parse options pertaining to total force calculation - virtual int init_total_force_params(std::string const &conf); - /// \brief After construction, set data related to dependency handling int setup(); - /// \brief Default constructor (used when \link colvar::cvc \endlink - /// objects are declared within other ones) - cvc(); - - /// Destructor - virtual ~cvc(); - - /// \brief Implementation of the feature list for colvar - static std::vector cvc_features; - /// \brief Implementation of the feature list accessor for colvar virtual const std::vector &features() const { return cvc_features; } + virtual std::vector &modify_features() { return cvc_features; } + static void delete_features() { for (size_t i=0; i < cvc_features.size(); i++) { delete cvc_features[i]; @@ -206,16 +177,12 @@ class colvar::cvc /// collective variable force, usually coming from the biases and /// eventually manipulated by the parent \link colvar \endlink /// object - virtual void apply_force(colvarvalue const &cvforce) = 0; + virtual void apply_force(colvarvalue const &cvforce); - /// \brief Square distance between x1 and x2 (can be redefined to - /// transparently implement constraints, symmetries and - /// periodicities) - /// - /// colvar::cvc::dist2() and the related functions are - /// declared as "const" functions, but not "static", because - /// additional parameters defining the metrics (e.g. the - /// periodicity) may be specific to each colvar::cvc object. + /// Square distance between x1 and x2 (can be redefined to transparently + /// implement metrics in multi-dimensional spaces with or without + /// constraints, symmetries and periodicities). The default implementation + /// assumes scalar numbers and no symmetries or periodicities. /// /// If symmetries or periodicities are present, the /// colvar::cvc::dist2() should be redefined to return the @@ -230,12 +197,6 @@ class colvar::cvc /// i.e. already deprived of its component normal to the constraint /// hypersurface. /// - /// Finally, another useful application, if you are performing very - /// many operations with these functions, could be to override the - /// \link colvarvalue \endlink member functions and access directly - /// its member data. For instance: to define dist2(x1,x2) as - /// (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in - /// case of a scalar \link colvarvalue \endlink type. virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; @@ -268,7 +229,7 @@ class colvar::cvc virtual int set_param(std::string const ¶m_name, void const *new_value); /// \brief Whether or not this CVC will be computed in parallel whenever possible - bool b_try_scalable; + bool b_try_scalable = true; /// Forcibly set value of CVC - useful for driving an external coordinate, /// eg. lambda dynamics @@ -278,9 +239,22 @@ class colvar::cvc protected: + /// Set the value of \link function_type \endlink and its dependencies + int set_function_type(std::string const &type); + /// Update the description string based on name and type int update_description(); + /// Parse a group definition + cvm::atom_group *parse_group(std::string const &conf, char const *group_key, + bool optional = false); + + /// \brief Parse options pertaining to total force calculation + virtual int init_total_force_params(std::string const &conf); + + /// \brief Implementation of the feature list for colvar + static std::vector cvc_features; + /// Record the type of this class as well as those it is derived from std::vector function_types; @@ -318,8 +292,8 @@ class colvar::cvc /// \brief Location of the upper boundary (not defined by user choice) colvarvalue upper_boundary; - /// \brief CVC-specific default colvar width - cvm::real width; + /// \brief CVC-specific default colvar width (default: not provided) + cvm::real width = 0.0; }; @@ -350,26 +324,19 @@ class colvar::distance { protected: /// First atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Second atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// Vector distance, cached to be recycled cvm::rvector dist_v; public: - distance(std::string const &conf); distance(); virtual ~distance() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -381,25 +348,22 @@ class colvar::distance_vec : public colvar::distance { public: - distance_vec(std::string const &conf); distance_vec(); virtual ~distance_vec() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - /// Redefined to handle the box periodicity - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the box periodicity - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the box periodicity - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; - /// \brief Colvar component: distance unit vector (direction) between /// centers of mass of two groups (colvarvalue::type_unit3vector type, /// range [-1:1]x[-1:1]x[-1:1]) @@ -407,21 +371,19 @@ class colvar::distance_dir : public colvar::distance { public: - distance_dir(std::string const &conf); distance_dir(); virtual ~distance_dir() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - /// Redefined to override the distance ones - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to override the distance ones - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to override the distance ones - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -433,11 +395,11 @@ class colvar::distance_z { protected: /// Main atom group - cvm::atom_group *main; + cvm::atom_group *main = nullptr; /// Reference atom group - cvm::atom_group *ref1; + cvm::atom_group *ref1 = nullptr; /// Optional, second ref atom group - cvm::atom_group *ref2; + cvm::atom_group *ref2 = nullptr; /// Vector on which the distance vector is projected cvm::rvector axis; /// Norm of the axis @@ -445,24 +407,15 @@ class colvar::distance_z /// Vector distance, cached to be recycled cvm::rvector dist_v; /// Flag: using a fixed axis vector? - bool fixed_axis; + bool fixed_axis = true; public: - distance_z(std::string const &conf); distance_z(); virtual ~distance_z() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// \brief Redefined to make use of the user-provided period - virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -478,20 +431,12 @@ class colvar::distance_xy /// Vector distances cvm::rvector v12, v13; public: - distance_xy(std::string const &conf); distance_xy(); virtual ~distance_xy() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -500,28 +445,16 @@ class colvar::distance_xy class colvar::polar_phi : public colvar::cvc { -public: - polar_phi(std::string const &conf); - polar_phi(); - virtual ~polar_phi() {} protected: - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; cvm::real r, theta, phi; + public: + polar_phi(); + virtual ~polar_phi() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - /// Redefined to handle the 2*PI periodicity - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -531,27 +464,18 @@ class colvar::polar_theta : public colvar::cvc { public: - polar_theta(std::string const &conf); polar_theta(); virtual ~polar_theta() {} + virtual int init(std::string const &conf); protected: - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; cvm::real r, theta, phi; public: virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - /// Redefined to override the distance ones - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to override the distance ones - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to override the distance ones - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; + /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power, /// as in NMR refinements(colvarvalue::type_scalar type, range (0:*)) class colvar::distance_inv @@ -559,23 +483,17 @@ class colvar::distance_inv { protected: /// First atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Second atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// Components of the distance vector orthogonal to the axis - int exponent; + int exponent = 6; public: - distance_inv(std::string const &conf); + distance_inv(); virtual ~distance_inv() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -587,45 +505,41 @@ class colvar::distance_pairs { protected: /// First atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Second atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; public: - distance_pairs(std::string const &conf); distance_pairs(); virtual ~distance_pairs() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); + /// Redefined to deal with multiple dimensions + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; - /// \brief Colvar component: dipole magnitude of a molecule class colvar::dipole_magnitude : public colvar::cvc { protected: /// Dipole atom group - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; cvm::atom_pos dipoleV; public: - /// Initialize by parsing the configuration - dipole_magnitude (std::string const &conf); - dipole_magnitude (cvm::atom const &a1); dipole_magnitude(); virtual ~dipole_magnitude() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - //virtual void calc_force_invgrads(); - //virtual void calc_Jacobian_derivative(); - virtual void apply_force (colvarvalue const &force); - virtual cvm::real dist2 (colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad (colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad (colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -637,21 +551,15 @@ class colvar::gyration { protected: /// Atoms involved - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; public: - gyration(std::string const &conf); + gyration(); virtual ~gyration() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -662,19 +570,10 @@ class colvar::inertia : public colvar::gyration { public: - /// Constructor - inertia(std::string const &conf); inertia(); virtual ~inertia() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -688,19 +587,11 @@ class colvar::inertia_z /// Vector on which the inertia tensor is projected cvm::rvector axis; public: - /// Constructor - inertia_z(std::string const &conf); inertia_z(); virtual ~inertia_z() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -713,7 +604,7 @@ class colvar::eigenvector protected: /// Atom group - cvm::atom_group * atoms; + cvm::atom_group * atoms = nullptr; /// Reference coordinates std::vector ref_pos; @@ -726,20 +617,13 @@ class colvar::eigenvector public: - /// Constructor - eigenvector(std::string const &conf); + eigenvector(); virtual ~eigenvector() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -752,11 +636,11 @@ class colvar::angle protected: /// Atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// Atom group - cvm::atom_group *group3; + cvm::atom_group *group3 = nullptr; /// Inter site vectors cvm::rvector r21, r23; @@ -768,25 +652,18 @@ class colvar::angle /// Compute total force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) /// (or to allow dummy atoms) - bool b_1site_force; + bool b_1site_force = false; public: - /// Initialize by parsing the configuration - angle(std::string const &conf); + angle(); /// \brief Initialize the three groups after three atoms angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3); virtual ~angle() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -799,11 +676,11 @@ class colvar::dipole_angle protected: /// Dipole atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// Atom group - cvm::atom_group *group3; + cvm::atom_group *group3 = nullptr; /// Inter site vectors cvm::rvector r21, r23; @@ -815,24 +692,14 @@ class colvar::dipole_angle /// Compute total force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) /// (or to allow dummy atoms) - bool b_1site_force; + bool b_1site_force = false; public: - /// Initialize by parsing the configuration - dipole_angle (std::string const &conf); - /// \brief Initialize the three groups after three atoms - dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3); dipole_angle(); virtual ~dipole_angle() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force (colvarvalue const &force); - virtual cvm::real dist2 (colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad (colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad (colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -845,45 +712,31 @@ class colvar::dihedral protected: /// Atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// Atom group - cvm::atom_group *group3; + cvm::atom_group *group3 = nullptr; /// Atom group - cvm::atom_group *group4; + cvm::atom_group *group4 = nullptr; /// Inter site vectors cvm::rvector r12, r23, r34; /// \brief Compute total force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) - bool b_1site_force; + bool b_1site_force = false; public: - /// Initialize by parsing the configuration - dihedral(std::string const &conf); /// \brief Initialize the four groups after four atoms dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4); dihedral(); virtual ~dihedral() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - - /// Redefined to handle the 2*PI periodicity - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -895,46 +748,39 @@ class colvar::coordnum { protected: /// First atom group - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// Second atom group - cvm::atom_group *group2; + cvm::atom_group *group2 = nullptr; /// \brief "Cutoff" for isotropic calculation (default) cvm::real r0; /// \brief "Cutoff vector" for anisotropic calculation cvm::rvector r0_vec; /// \brief Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used - bool b_anisotropic; + bool b_anisotropic = false; /// Integer exponent of the function numerator - int en; + int en = 6; /// Integer exponent of the function denominator - int ed; + int ed = 12; /// If true, group2 will be treated as a single atom - bool b_group2_center_only; + bool b_group2_center_only = false; /// Tolerance for the pair list - cvm::real tolerance; + cvm::real tolerance = 0.0; /// Frequency of update of the pair list - int pairlist_freq; + int pairlist_freq = 100; /// Pair list - bool *pairlist; + bool *pairlist = nullptr; public: - coordnum(std::string const &conf); - ~coordnum(); - + coordnum(); + virtual ~coordnum(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; enum { ef_null = 0, @@ -979,31 +825,25 @@ class colvar::selfcoordnum protected: /// Selected atoms - cvm::atom_group *group1; + cvm::atom_group *group1 = nullptr; /// \brief "Cutoff" for isotropic calculation (default) cvm::real r0; /// Integer exponent of the function numerator - int en; + int en = 6; /// Integer exponent of the function denominator - int ed; - cvm::real tolerance; - int pairlist_freq; - bool *pairlist; + int ed = 12; + cvm::real tolerance = 0.0; + int pairlist_freq = 100; + + bool *pairlist = nullptr; public: - selfcoordnum(std::string const &conf); + selfcoordnum(); ~selfcoordnum(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; /// Main workhorse function template int compute_selfcoordnum(); @@ -1023,25 +863,17 @@ class colvar::groupcoordnum cvm::rvector r0_vec; /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be /// used - bool b_anisotropic; + bool b_anisotropic = false; /// Integer exponent of the function numerator - int en; + int en = 6; /// Integer exponent of the function denominator - int ed; + int ed = 12; public: - /// Constructor - groupcoordnum(std::string const &conf); + groupcoordnum(); virtual ~groupcoordnum() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1056,73 +888,22 @@ class colvar::h_bond /// \brief "Cutoff" distance between acceptor and donor cvm::real r0; /// Integer exponent of the function numerator - int en; + int en = 6; /// Integer exponent of the function denominator - int ed; + int ed = 8; public: - h_bond(std::string const &conf); /// Constructor for atoms already allocated h_bond(cvm::atom const &acceptor, cvm::atom const &donor, cvm::real r0, int en, int ed); h_bond(); virtual ~h_bond() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; - -/// \brief Colvar component: alpha helix content of a contiguous -/// segment of 5 or more residues, implemented as a sum of phi/psi -/// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type, -/// range [0:1]) -// class colvar::alpha_dihedrals -// : public colvar::cvc -// { -// protected: - -// /// Alpha-helical reference phi value -// cvm::real phi_ref; - -// /// Alpha-helical reference psi value -// cvm::real psi_ref; - -// /// List of phi dihedral angles -// std::vector phi; - -// /// List of psi dihedral angles -// std::vector psi; - -// /// List of hydrogen bonds -// std::vector hb; - -// public: - -// alpha_dihedrals (std::string const &conf); -// alpha_dihedrals(); -// virtual ~alpha_dihedrals() {} -// virtual void calc_value(); -// virtual void calc_gradients(); -// virtual void apply_force (colvarvalue const &force); -// virtual cvm::real dist2 (colvarvalue const &x1, -// colvarvalue const &x2) const; -// virtual colvarvalue dist2_lgrad (colvarvalue const &x1, -// colvarvalue const &x2) const; -// virtual colvarvalue dist2_rgrad (colvarvalue const &x1, -// colvarvalue const &x2) const; -// }; - - - /// \brief Colvar component: alpha helix content of a contiguous /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca /// angles and hydrogen bonds (colvarvalue::type_scalar type, range @@ -1133,10 +914,10 @@ class colvar::alpha_angles protected: /// Reference Calpha-Calpha angle (default: 88 degrees) - cvm::real theta_ref; + cvm::real theta_ref = 88.0; /// Tolerance on the Calpha-Calpha angle - cvm::real theta_tol; + cvm::real theta_tol = 15.0; /// List of Calpha-Calpha angles std::vector theta; @@ -1144,25 +925,28 @@ class colvar::alpha_angles /// List of hydrogen bonds std::vector hb; - /// Contribution of the hb terms - cvm::real hb_coeff; + /// Contribution of the HB terms + cvm::real hb_coeff = 0.5; + + /// Cutoff for HB + cvm::real r0; + + /// Integer exponent of the HB numerator + int en = 6; + + /// Integer exponent of the HB denominator + int ed = 8; public: - alpha_angles(std::string const &conf); alpha_angles(); virtual ~alpha_angles(); + virtual int init(std::string const &conf); void calc_value(); void calc_gradients(); /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1181,20 +965,14 @@ class colvar::dihedPC public: - dihedPC(std::string const &conf); dihedPC(); - virtual ~dihedPC(); - void calc_value(); - void calc_gradients(); + virtual ~dihedPC(); + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs - void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); - void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; + virtual void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); + virtual void apply_force(colvarvalue const &force); }; @@ -1209,7 +987,7 @@ class colvar::orientation protected: /// Atom group - cvm::atom_group * atoms; + cvm::atom_group * atoms = nullptr; /// Center of geometry of the group cvm::atom_pos atoms_cog; @@ -1232,43 +1010,44 @@ class colvar::orientation public: - orientation(std::string const &conf); orientation(); - virtual int init(std::string const &conf); virtual ~orientation(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; + /// Redefined to use quaternion metrics + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use quaternion metrics + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use quaternion metrics + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use quaternion metrics + virtual void wrap(colvarvalue &x_unwrapped) const; }; - -/// \brief Colvar component: angle of rotation with respect to a set -/// of reference coordinates (colvarvalue::type_scalar type, range -/// [0:PI)) +/// Colvar component: angle of rotation with respect to a set of reference coordinates +/// (colvarvalue::type_scalar type, range [0:PI)) +/// This is also used to derive all other sub-rotation variables that return a scalar value class colvar::orientation_angle : public colvar::orientation { public: - orientation_angle(std::string const &conf); - virtual int init(std::string const &conf); + orientation_angle(); virtual ~orientation_angle() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; + /// Redefined to use scalar metrics + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use scalar metrics + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use scalar metrics + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use scalar metrics + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1277,23 +1056,14 @@ class colvar::orientation_angle /// of reference coordinates (colvarvalue::type_scalar type, range /// [-1:1]) class colvar::orientation_proj - : public colvar::orientation + : public colvar::orientation_angle { public: - orientation_proj(std::string const &conf); orientation_proj(); - virtual int init(std::string const &conf); virtual ~orientation_proj() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1301,7 +1071,7 @@ class colvar::orientation_proj /// \brief Colvar component: projection of the orientation vector onto /// a predefined axis (colvarvalue::type_scalar type, range [-1:1]) class colvar::tilt - : public colvar::orientation + : public colvar::orientation_proj { protected: @@ -1309,18 +1079,11 @@ class colvar::tilt public: - tilt(std::string const &conf); - virtual int init(std::string const &conf); + tilt(); virtual ~tilt() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1328,98 +1091,47 @@ class colvar::tilt /// \brief Colvar component: angle of rotation around a predefined /// axis (colvarvalue::type_scalar type, range [-PI:PI]) class colvar::spin_angle - : public colvar::orientation + : public colvar::tilt { -protected: - - cvm::rvector axis; - public: - spin_angle(std::string const &conf); spin_angle(); - virtual int init(std::string const &conf); virtual ~spin_angle() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - /// Redefined to handle the 2*PI periodicity - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x_unwrapped) const; }; class colvar::euler_phi - : public colvar::orientation + : public colvar::orientation_angle { public: - euler_phi(std::string const &conf); euler_phi(); - virtual int init(std::string const &conf); virtual ~euler_phi() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x_unwrapped) const; }; class colvar::euler_psi - : public colvar::orientation + : public colvar::orientation_angle { public: - euler_psi(std::string const &conf); euler_psi(); - virtual int init(std::string const &conf); virtual ~euler_psi() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x_unwrapped) const; }; class colvar::euler_theta - : public colvar::orientation + : public colvar::orientation_angle { public: - euler_theta(std::string const &conf); euler_theta(); - virtual int init(std::string const &conf); virtual ~euler_theta() {} virtual void calc_value(); virtual void calc_gradients(); - virtual void apply_force(colvarvalue const &force); - // theta angle is a scalar variable and not periodic - // we need to override the virtual functions from orientation - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1433,33 +1145,29 @@ class colvar::rmsd protected: /// Atom group - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; /// Reference coordinates (for RMSD calculation only) /// Includes sets with symmetry permutations (n_permutations * n_atoms) std::vector ref_pos; /// Number of permutations of symmetry-related atoms - size_t n_permutations; + size_t n_permutations = 1; /// Index of the permutation yielding the smallest RMSD (0 for identity) - size_t best_perm_index; -public: + size_t best_perm_index = 0; - /// Constructor - rmsd(std::string const &conf); + /// Permutation RMSD input parsing + int init_permutation(std::string const &conf); + +public: + rmsd(); virtual ~rmsd() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); - virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1471,16 +1179,27 @@ class colvar::cartesian { protected: /// Atom group - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; /// Which Cartesian coordinates to include std::vector axes; public: - cartesian(std::string const &conf); cartesian(); virtual ~cartesian() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); + /// Redefined to deal with multiple dimensions + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to deal with multiple dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1492,18 +1211,11 @@ class colvar::alch_lambda protected: // No atom groups needed public: - alch_lambda(std::string const &conf); alch_lambda(); virtual ~alch_lambda() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; @@ -1515,36 +1227,14 @@ class colvar::alch_Flambda protected: // No atom groups needed public: - alch_Flambda(std::string const &conf); alch_Flambda(); virtual ~alch_Flambda() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); - virtual cvm::real dist2(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const; - virtual colvarvalue dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const; }; -class colvar::componentDisabled - : public colvar::cvc -{ -public: - componentDisabled(std::string const & /* conf */) { - cvm::error("Error: this component is not enabled in the current build; please see https://colvars.github.io/README-c++11.html"); - } - virtual ~componentDisabled() {} - virtual void calc_value() {} - virtual void calc_gradients() {} - virtual void apply_force(colvarvalue const & /* force */) {} -}; - - - class colvar::CartesianBasedPath : public colvar::cvc { @@ -1552,21 +1242,23 @@ class colvar::CartesianBasedPath virtual void computeDistanceBetweenReferenceFrames(std::vector& result); virtual void computeDistanceToReferenceFrames(std::vector& result); /// Selected atoms - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; /// Fitting options - bool has_user_defined_fitting; + bool has_user_defined_fitting = false; /// Reference frames std::vector> reference_frames; std::vector> reference_fitting_frames; /// Atom groups for RMSD calculation together with reference frames std::vector comp_atoms; /// Total number of reference frames - size_t total_reference_frames; + size_t total_reference_frames = 0; public: - CartesianBasedPath(std::string const &conf); + CartesianBasedPath(); virtual ~CartesianBasedPath(); + virtual int init(std::string const &conf); virtual void calc_value() = 0; - virtual void apply_force(colvarvalue const &force) = 0; + /// Redefined to raise error because this is an abstract type + virtual void apply_force(colvarvalue const &force); }; /// \brief Colvar component: alternative path collective variable using geometry, variable s @@ -1582,8 +1274,9 @@ class colvar::gspath virtual void prepareVectors(); virtual void updateDistanceToReferenceFrames(); public: - gspath(std::string const &conf); + gspath(); virtual ~gspath() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1604,8 +1297,9 @@ class colvar::gzpath virtual void prepareVectors(); virtual void updateDistanceToReferenceFrames(); public: - gzpath(std::string const &conf); + gzpath(); virtual ~gzpath() {} + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1623,11 +1317,23 @@ class colvar::linearCombination protected: cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const; public: - linearCombination(std::string const &conf); + linearCombination(); virtual ~linearCombination(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); + /// Redefined to allow arbitrary dimensions + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1636,7 +1342,7 @@ class colvar::customColvar : public colvar::linearCombination { protected: - bool use_custom_function; + bool use_custom_function = false; #ifdef LEPTON /// Vector of evaluators for custom functions using Lepton std::vector value_evaluators; @@ -1646,11 +1352,12 @@ class colvar::customColvar std::vector value_eval_var_refs; std::vector grad_eval_var_refs; /// Unused value that is written to when a variable simplifies out of a Lepton expression - double dev_null; + double dev_null = 0.0; #endif public: - customColvar(std::string const &conf); + customColvar(); virtual ~customColvar(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1668,17 +1375,27 @@ class colvar::CVBasedPath /// If all sub-cvs use explicit gradients then we also use it bool use_explicit_gradients; /// Total number of reference frames - size_t total_reference_frames; + size_t total_reference_frames = 0; protected: virtual void computeDistanceToReferenceFrames(std::vector& result); /// Helper function to determine the distance between reference frames virtual void computeDistanceBetweenReferenceFrames(std::vector& result) const; cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const; public: - CVBasedPath(std::string const &conf); + CVBasedPath(); virtual ~CVBasedPath(); + virtual int init(std::string const &conf); virtual void calc_value() = 0; - virtual void apply_force(colvarvalue const &force) = 0; + /// Redefined to raise error because this is an abstract type + virtual void apply_force(colvarvalue const &force); + /// Redefined to use the metric of the returned colvarvalue (defined at runtime) + virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use the metric of the returned colvarvalue (defined at runtime) + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use the metric of the returned colvarvalue (defined at runtime) + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; + /// Redefined to use the metric of the returned colvarvalue (defined at runtime) + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1693,8 +1410,9 @@ class colvar::gspathCV virtual void updateDistanceToReferenceFrames(); virtual void prepareVectors(); public: - gspathCV(std::string const &conf); + gspathCV(); virtual ~gspathCV(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1709,8 +1427,9 @@ class colvar::gzpathCV virtual void updateDistanceToReferenceFrames(); virtual void prepareVectors(); public: - gzpathCV(std::string const &conf); + gzpathCV(); virtual ~gzpathCV(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1725,8 +1444,9 @@ class colvar::aspath std::unique_ptr impl_; friend struct ArithmeticPathImpl; public: - aspath(std::string const &conf); + aspath(); virtual ~aspath(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1739,8 +1459,9 @@ class colvar::azpath std::unique_ptr impl_; friend struct ArithmeticPathImpl; public: - azpath(std::string const &conf); + azpath(); virtual ~azpath(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1753,8 +1474,9 @@ class colvar::aspathCV std::unique_ptr impl_; friend struct ArithmeticPathImpl; public: - aspathCV(std::string const &conf); + aspathCV(); virtual ~aspathCV(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1768,8 +1490,9 @@ class colvar::azpathCV std::unique_ptr impl_; friend struct ArithmeticPathImpl; public: - azpathCV(std::string const &conf); + azpathCV(); virtual ~azpathCV(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1788,13 +1511,25 @@ class colvar::neuralNetwork /// actual computation happens in neuralnetworkCV::neuralNetworkCompute std::unique_ptr nn; /// the index of nn output components - size_t m_output_index; + size_t m_output_index = 0; public: - neuralNetwork(std::string const &conf); + neuralNetwork(); virtual ~neuralNetwork(); + virtual int init(std::string const &conf); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); + /// Redefined to allow arbitrary dimensions + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to allow arbitrary dimensions + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1806,7 +1541,6 @@ class colvar::map_total public: map_total(); - map_total(std::string const &conf); virtual ~map_total() {} virtual int init(std::string const &conf); virtual void calc_value(); @@ -1819,13 +1553,13 @@ class colvar::map_total std::string volmap_name; /// Numeric identifier of the map object (as used by the simulation engine) - int volmap_id; + int volmap_id = -1; /// Index of the map objet in the proxy arrays - int volmap_index; + int volmap_index = -1; /// Group of atoms selected internally (optional) - cvm::atom_group *atoms; + cvm::atom_group *atoms = nullptr; /// Weights assigned to each atom (default: uniform weights) std::vector atom_weights; @@ -1833,36 +1567,4 @@ class colvar::map_total -// metrics functions for cvc implementations - -// simple definitions of the distance functions; these are useful only -// for optimization (the type check performed in the default -// colvarcomp functions is skipped) - -// definitions assuming the scalar type - -#define simple_scalar_dist_functions(TYPE) \ - \ - \ - cvm::real colvar::TYPE::dist2(colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \ - } \ - \ - \ - colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return 2.0 * (x1.real_value - x2.real_value); \ - } \ - \ - \ - colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return this->dist2_lgrad(x2, x1); \ - } \ - - #endif diff --git a/src/colvarcomp_alchlambda.cpp b/src/colvarcomp_alchlambda.cpp index dbb1ece45..a175ea43e 100644 --- a/src/colvarcomp_alchlambda.cpp +++ b/src/colvarcomp_alchlambda.cpp @@ -16,8 +16,7 @@ #include "colvarcomp.h" -colvar::alch_lambda::alch_lambda(std::string const &conf) - : cvc(conf) +colvar::alch_lambda::alch_lambda() { set_function_type("alchLambda"); @@ -55,12 +54,9 @@ void colvar::alch_lambda::apply_force(colvarvalue const & /* force */) cvm::proxy->set_alch_lambda(x.real_value); } -simple_scalar_dist_functions(alch_lambda) - -colvar::alch_Flambda::alch_Flambda(std::string const &conf) - : cvc(conf) +colvar::alch_Flambda::alch_Flambda() { set_function_type("alch_Flambda"); @@ -102,4 +98,3 @@ void colvar::alch_Flambda::apply_force(colvarvalue const &force) cvm::proxy->indirect_lambda_biasing_force += d2E_dlambda2 * f; } -simple_scalar_dist_functions(alch_Flambda) diff --git a/src/colvarcomp_angles.cpp b/src/colvarcomp_angles.cpp index 6e02f2c9c..56894e9f5 100644 --- a/src/colvarcomp_angles.cpp +++ b/src/colvarcomp_angles.cpp @@ -12,8 +12,7 @@ #include "colvarcomp.h" -colvar::angle::angle(std::string const &conf) - : cvc(conf) +colvar::angle::angle() { set_function_type("angle"); init_as_angle(); @@ -21,26 +20,25 @@ colvar::angle::angle(std::string const &conf) provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); +} + + +int colvar::angle::init(std::string const &conf) +{ + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); group3 = parse_group(conf, "group3"); - init_total_force_params(conf); + error_code |= init_total_force_params(conf); + + return error_code; } -colvar::angle::angle(cvm::atom const &a1, - cvm::atom const &a2, - cvm::atom const &a3) +colvar::angle::angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3) : angle() { - set_function_type("angle"); - init_as_angle(); - - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); - group1 = new cvm::atom_group(std::vector(1, a1)); group2 = new cvm::atom_group(std::vector(1, a2)); group3 = new cvm::atom_group(std::vector(1, a3)); @@ -120,57 +118,24 @@ void colvar::angle::calc_Jacobian_derivative() } -void colvar::angle::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); - - if (!group3->noforce) - group3->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(angle) - - -colvar::dipole_angle::dipole_angle(std::string const &conf) - : cvc(conf) +colvar::dipole_angle::dipole_angle() { set_function_type("dipoleAngle"); init_as_angle(); - - group1 = parse_group(conf, "group1"); - group2 = parse_group(conf, "group2"); - group3 = parse_group(conf, "group3"); - - init_total_force_params(conf); } -colvar::dipole_angle::dipole_angle(cvm::atom const &a1, - cvm::atom const &a2, - cvm::atom const &a3) +int colvar::dipole_angle::init(std::string const &conf) { - set_function_type("dipoleAngle"); - init_as_angle(); - - group1 = new cvm::atom_group(std::vector(1, a1)); - group2 = new cvm::atom_group(std::vector(1, a2)); - group3 = new cvm::atom_group(std::vector(1, a3)); - register_atom_group(group1); - register_atom_group(group2); - register_atom_group(group3); -} + int error_code = cvc::init(conf); + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); + group3 = parse_group(conf, "group3"); -colvar::dipole_angle::dipole_angle() -{ - set_function_type("dipoleAngle"); - init_as_angle(); + error_code |= init_total_force_params(conf); + return error_code; } @@ -229,52 +194,35 @@ void colvar::dipole_angle::calc_gradients() } -void colvar::dipole_angle::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); - - if (!group3->noforce) - group3->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(dipole_angle) - - -colvar::dihedral::dihedral(std::string const &conf) - : cvc(conf) +colvar::dihedral::dihedral() { set_function_type("dihedral"); init_as_periodic_angle(); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); +} + + +int colvar::dihedral::init(std::string const &conf) +{ + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); group3 = parse_group(conf, "group3"); group4 = parse_group(conf, "group4"); - init_total_force_params(conf); + error_code |= init_total_force_params(conf); + return error_code; } -colvar::dihedral::dihedral(cvm::atom const &a1, - cvm::atom const &a2, - cvm::atom const &a3, +colvar::dihedral::dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4) + : dihedral() { - set_function_type("dihedral"); - init_as_periodic_angle(); - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); - b_1site_force = false; group1 = new cvm::atom_group(std::vector(1, a1)); @@ -288,16 +236,6 @@ colvar::dihedral::dihedral(cvm::atom const &a1, } -colvar::dihedral::dihedral() -{ - set_function_type("dihedral"); - init_as_periodic_angle(); - enable(f_cvc_periodic); - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); -} - - void colvar::dihedral::calc_value() { cvm::atom_pos const g1_pos = group1->center_of_mass(); @@ -323,7 +261,7 @@ void colvar::dihedral::calc_value() cvm::real const sin_phi = n1 * r34 * r23.norm(); x.real_value = (180.0/PI) * cvm::atan2(sin_phi, cos_phi); - this->wrap(x); + wrap(x); } @@ -439,81 +377,22 @@ void colvar::dihedral::calc_Jacobian_derivative() } -void colvar::dihedral::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); - - if (!group3->noforce) - group3->apply_colvar_force(force.real_value); - - if (!group4->noforce) - group4->apply_colvar_force(force.real_value); -} - - -// metrics functions for cvc implementations with a periodicity - -cvm::real colvar::dihedral::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return diff * diff; -} - - -colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return 2.0 * diff; -} - - -colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return (-2.0) * diff; -} - - -void colvar::dihedral::wrap(colvarvalue &x_unwrapped) const -{ - if ((x_unwrapped.real_value - wrap_center) >= 180.0) { - x_unwrapped.real_value -= 360.0; - return; - } - if ((x_unwrapped.real_value - wrap_center) < -180.0) { - x_unwrapped.real_value += 360.0; - return; - } -} - - -colvar::polar_theta::polar_theta(std::string const &conf) - : cvc(conf) +colvar::polar_theta::polar_theta() { + r = theta = phi = 0.0; set_function_type("polarTheta"); enable(f_cvc_com_based); - - atoms = parse_group(conf, "atoms"); - init_total_force_params(conf); - x.type(colvarvalue::type_scalar); + init_as_angle(); } -colvar::polar_theta::polar_theta() +int colvar::polar_theta::init(std::string const &conf) { - set_function_type("polarTheta"); - x.type(colvarvalue::type_scalar); + int error_code = cvc::init(conf); + atoms = parse_group(conf, "atoms"); + error_code |= init_total_force_params(conf); + return error_code; } @@ -540,32 +419,22 @@ void colvar::polar_theta::calc_gradients() } -void colvar::polar_theta::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(polar_theta) - -colvar::polar_phi::polar_phi(std::string const &conf) - : cvc(conf) +colvar::polar_phi::polar_phi() { + r = theta = phi = 0.0; set_function_type("polarPhi"); - init_as_periodic_angle(); enable(f_cvc_com_based); - - atoms = parse_group(conf, "atoms"); - init_total_force_params(conf); + init_as_periodic_angle(); } -colvar::polar_phi::polar_phi() +int colvar::polar_phi::init(std::string const &conf) { - set_function_type("polarPhi"); - init_as_periodic_angle(); + int error_code = cvc::init(conf); + atoms = parse_group(conf, "atoms"); + error_code |= init_total_force_params(conf); + return error_code; } @@ -587,55 +456,3 @@ void colvar::polar_phi::calc_gradients() (180.0/PI) * cvm::cos(phi) / (r*cvm::sin(theta)), 0.)); } - - -void colvar::polar_phi::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - -// Same as dihedral, for polar_phi - -cvm::real colvar::polar_phi::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return diff * diff; -} - - -colvarvalue colvar::polar_phi::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return 2.0 * diff; -} - - -colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return (-2.0) * diff; -} - - -void colvar::polar_phi::wrap(colvarvalue &x_unwrapped) const -{ - if ((x_unwrapped.real_value - wrap_center) >= 180.0) { - x_unwrapped.real_value -= 360.0; - return; - } - - if ((x_unwrapped.real_value - wrap_center) < -180.0) { - x_unwrapped.real_value += 360.0; - return; - } - - return; -} diff --git a/src/colvarcomp_apath.cpp b/src/colvarcomp_apath.cpp index 1a8b94550..9e8256957 100644 --- a/src/colvarcomp_apath.cpp +++ b/src/colvarcomp_apath.cpp @@ -134,17 +134,27 @@ struct ArithmeticPathImpl: public ArithmeticPathCV::ArithmeticPathBasesize(); std::vector p_weights(num_atoms, std::sqrt(1.0 / num_atoms)); // ArithmeticPathCV::ArithmeticPathBase::initialize(num_atoms, total_reference_frames, p_lambda, reference_frames[0], p_weights); + if (impl_) impl_.reset(); impl_ = std::unique_ptr(new ArithmeticPathImpl(num_atoms, total_reference_frames, p_lambda, p_weights)); cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n")); + return error_code; } colvar::aspath::~aspath() {} @@ -180,16 +190,26 @@ void colvar::aspath::apply_force(colvarvalue const &force) { } } -colvar::azpath::azpath(std::string const &conf): CartesianBasedPath(conf) { - function_type = "azpath"; +colvar::azpath::azpath() +{ + set_function_type("azpath"); + x.type(colvarvalue::type_scalar); +} + +int colvar::azpath::init(std::string const &conf) +{ + int error_code = CartesianBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); x.type(colvarvalue::type_scalar); cvm::real p_lambda; get_keyval(conf, "lambda", p_lambda, -1.0); const size_t num_atoms = atoms->size(); std::vector p_weights(num_atoms, std::sqrt(1.0 / num_atoms)); + if (impl_) impl_.reset(); impl_ = std::unique_ptr(new ArithmeticPathImpl(num_atoms, total_reference_frames, p_lambda, p_weights)); cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n")); + return error_code; } colvar::azpath::~azpath() {} @@ -225,15 +245,23 @@ void colvar::azpath::apply_force(colvarvalue const &force) { } } -colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) { +colvar::aspathCV::aspathCV() +{ set_function_type("aspathCV"); + x.type(colvarvalue::type_scalar); +} + +int colvar::aspathCV::init(std::string const &conf) +{ + int error_code = CVBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); std::vector p_weights(cv.size(), 1.0); get_keyval(conf, "weights", p_weights, std::vector(cv.size(), 1.0)); - x.type(colvarvalue::type_scalar); use_explicit_gradients = true; cvm::real p_lambda; get_keyval(conf, "lambda", p_lambda, -1.0); + if (impl_) impl_.reset(); impl_ = std::unique_ptr(new ArithmeticPathImpl(cv.size(), total_reference_frames, p_lambda, p_weights)); cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n")); for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { @@ -242,6 +270,7 @@ colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) { } cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(p_weights[i_cv]) + std::string("\n")); } + return error_code; } colvar::aspathCV::~aspathCV() {} @@ -321,15 +350,23 @@ void colvar::aspathCV::apply_force(colvarvalue const &force) { } } -colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) { +colvar::azpathCV::azpathCV() +{ set_function_type("azpathCV"); + x.type(colvarvalue::type_scalar); +} + +int colvar::azpathCV::init(std::string const &conf) +{ + int error_code = CVBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); std::vector p_weights(cv.size(), 1.0); get_keyval(conf, "weights", p_weights, std::vector(cv.size(), 1.0)); - x.type(colvarvalue::type_scalar); use_explicit_gradients = true; cvm::real p_lambda; get_keyval(conf, "lambda", p_lambda, -1.0); + if (impl_) impl_.reset(); impl_ = std::unique_ptr(new ArithmeticPathImpl(cv.size(), total_reference_frames, p_lambda, p_weights)); cvm::log(std::string("Lambda is ") + cvm::to_str(impl_->get_lambda()) + std::string("\n")); for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { @@ -338,6 +375,7 @@ colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) { } cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(p_weights[i_cv]) + std::string("\n")); } + return error_code; } void colvar::azpathCV::calc_value() { diff --git a/src/colvarcomp_combination.cpp b/src/colvarcomp_combination.cpp index 974cc870c..01c7ffc32 100644 --- a/src/colvarcomp_combination.cpp +++ b/src/colvarcomp_combination.cpp @@ -9,14 +9,26 @@ #include "colvarcomp.h" -colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) { + +colvar::linearCombination::linearCombination() +{ + set_function_type("linearCombination"); +} + + +int colvar::linearCombination::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + if (error_code != COLVARS_OK) return error_code; + // Lookup all available sub-cvcs for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) { if (key_lookup(conf, it_cv_map->first.c_str())) { std::vector sub_cvc_confs; get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs); for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) { - cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf))); + cv.push_back((it_cv_map->second)()); + cv.back()->init(*(it_sub_cvc_conf)); } } } @@ -29,9 +41,8 @@ colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) } // Show useful error messages and prevent crashes if no sub CVC is found if (cv.size() == 0) { - cvm::error("Error: the CV " + name + - " expects one or more nesting components.\n"); - return; + return cvm::error("Error: the CV " + name + " expects one or more nesting components.\n", + COLVARS_INPUT_ERROR); } else { x.type(cv[0]->value()); x.reset(); @@ -45,6 +56,7 @@ colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) if (!use_explicit_gradients) { disable(f_cvc_explicit_gradient); } + return error_code; } cvm::real colvar::linearCombination::getPolynomialFactorOfCVGradient(size_t i_cv) const { @@ -124,8 +136,42 @@ void colvar::linearCombination::apply_force(colvarvalue const &force) { } } -colvar::customColvar::customColvar(std::string const &conf): linearCombination(conf) { - use_custom_function = false; + +cvm::real colvar::linearCombination::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x1.dist2(x2); +} + + +colvarvalue colvar::linearCombination::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return x1.dist2_grad(x2); +} + + +colvarvalue colvar::linearCombination::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return x2.dist2_grad(x1); +} + + +void colvar::linearCombination::wrap(colvarvalue & /* x_unwrapped */) const {} + + + +colvar::customColvar::customColvar() +{ + set_function_type("customColvar"); +} + + +int colvar::customColvar::init(std::string const &conf) +{ + int error_code = linearCombination::init(conf); + if (error_code != COLVARS_OK) return error_code; + // code swipe from colvar::init_custom_function std::string expr_in, expr; size_t pos = 0; // current position in config string @@ -146,7 +192,7 @@ colvar::customColvar::customColvar(std::string const &conf): linearCombination(c pexpr = Lepton::Parser::parse(expr); pexprs.push_back(pexpr); } catch (...) { - cvm::error("Error parsing expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); + return cvm::error("Error parsing expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); } try { value_evaluators.push_back(new Lepton::CompiledExpression(pexpr.createCompiledExpression())); @@ -164,7 +210,7 @@ colvar::customColvar::customColvar(std::string const &conf): linearCombination(c } } } catch (...) { - cvm::error("Error compiling expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); + return cvm::error("Error compiling expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR); } } while (key_lookup(conf, "customFunction", &expr_in, &pos)); // Now define derivative with respect to each scalar sub-component @@ -189,7 +235,7 @@ colvar::customColvar::customColvar(std::string const &conf): linearCombination(c } } if (value_evaluators.size() == 0) { - cvm::error("Error: no custom function defined.\n", COLVARS_INPUT_ERROR); + return cvm::error("Error: no custom function defined.\n", COLVARS_INPUT_ERROR); } if (value_evaluators.size() != 1) { x.type(colvarvalue::type_vector); @@ -197,14 +243,17 @@ colvar::customColvar::customColvar(std::string const &conf): linearCombination(c x.type(colvarvalue::type_scalar); } #else - cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" - "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", - COLVARS_INPUT_ERROR); + return cvm::error( + "customFunction requires the Lepton library, but it is not enabled during compilation.\n" + "Please refer to the Compilation Notes section of the Colvars manual for more " + "information.\n", + COLVARS_NOT_IMPLEMENTED); #endif } else { cvm::log("Warning: no customFunction specified.\n"); cvm::log("Warning: use linear combination instead.\n"); } + return error_code; } colvar::customColvar::~customColvar() { diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp index 66a45594c..5104b85eb 100644 --- a/src/colvarcomp_coordnums.cpp +++ b/src/colvarcomp_coordnums.cpp @@ -89,48 +89,47 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0, } -colvar::coordnum::coordnum(std::string const &conf) - : cvc(conf), b_anisotropic(false), pairlist(NULL) - +colvar::coordnum::coordnum() { set_function_type("coordNum"); x.type(colvarvalue::type_scalar); - colvarproxy *proxy = cvm::main()->proxy; + r0 = proxy->angstrom_to_internal(4.0); + r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0)); +} + + +int colvar::coordnum::init(std::string const &conf) +{ + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); - if (group1 == NULL || group2 == NULL) { - cvm::error("Error: failed to initialize atom groups.\n", - COLVARS_INPUT_ERROR); - return; + if (!group1 || !group2) { + return cvm::error("Error: failed to initialize atom groups.\n", COLVARS_INPUT_ERROR); } if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) { - cvm::error("Error: group1 and group2 share a common atom (number: " + - cvm::to_str(atom_number) + ")\n", COLVARS_INPUT_ERROR); - return; + error_code |= cvm::error( + "Error: group1 and group2 share a common atom (number: " + cvm::to_str(atom_number) + ")\n", + COLVARS_INPUT_ERROR); } if (group1->b_dummy) { - cvm::error("Error: only group2 is allowed to be a dummy atom\n", - COLVARS_INPUT_ERROR); - return; + error_code |= + cvm::error("Error: only group2 is allowed to be a dummy atom\n", COLVARS_INPUT_ERROR); } - bool const b_isotropic = get_keyval(conf, "cutoff", r0, - cvm::real(proxy->angstrom_to_internal(4.0))); + bool const b_isotropic = get_keyval(conf, "cutoff", r0, r0); - if (get_keyval(conf, "cutoff3", r0_vec, - cvm::rvector(proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0)))) { + if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) { if (b_isotropic) { - cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" " - "at the same time.\n", - COLVARS_INPUT_ERROR); - return; + error_code |= cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" " + "at the same time.\n", + COLVARS_INPUT_ERROR); } b_anisotropic = true; @@ -140,17 +139,17 @@ colvar::coordnum::coordnum(std::string const &conf) if (r0_vec.z < 0.0) r0_vec.z *= -1.0; } - get_keyval(conf, "expNumer", en, 6); - get_keyval(conf, "expDenom", ed, 12); + get_keyval(conf, "expNumer", en, en); + get_keyval(conf, "expDenom", ed, ed); if ( (en%2) || (ed%2) ) { - cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", + COLVARS_INPUT_ERROR); } if ( (en <= 0) || (ed <= 0) ) { - cvm::error("Error: negative exponent(s) provided.\n", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: negative exponent(s) provided.\n", + COLVARS_INPUT_ERROR); } if (!is_enabled(f_cvc_pbc_minimum_image)) { @@ -159,14 +158,14 @@ colvar::coordnum::coordnum(std::string const &conf) get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy); - get_keyval(conf, "tolerance", tolerance, 0.0); + get_keyval(conf, "tolerance", tolerance, tolerance); if (tolerance > 0) { cvm::main()->cite_feature("coordNum pairlist"); - get_keyval(conf, "pairListFrequency", pairlist_freq, 100); + get_keyval(conf, "pairListFrequency", pairlist_freq, pairlist_freq); if ( ! (pairlist_freq > 0) ) { - cvm::error("Error: non-positive pairlistfrequency provided.\n", - COLVARS_INPUT_ERROR); - return; // and do not allocate the pairlists below + return cvm::error("Error: non-positive pairlistfrequency provided.\n", + COLVARS_INPUT_ERROR); + // return and do not allocate the pairlists below } if (b_group2_center_only) { pairlist = new bool[group1->size()]; @@ -180,12 +179,14 @@ colvar::coordnum::coordnum(std::string const &conf) static_cast(group1->size()) : static_cast(group1->size() * group2->size())); + + return error_code; } colvar::coordnum::~coordnum() { - if (pairlist != NULL) { + if (pairlist) { delete [] pairlist; } } @@ -284,25 +285,23 @@ void colvar::coordnum::calc_gradients() } -void colvar::coordnum::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(coordnum) +// h_bond member functions +colvar::h_bond::h_bond() +{ + colvarproxy *proxy = cvm::main()->proxy; + r0 = proxy->angstrom_to_internal(3.3); + set_function_type("hBond"); + x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 1.0); +} -// h_bond member functions -colvar::h_bond::h_bond(std::string const &conf) -: cvc(conf) +int colvar::h_bond::init(std::string const &conf) { + int error_code = cvc::init(conf); + if (cvm::debug()) cvm::log("Initializing h_bond object.\n"); @@ -310,15 +309,12 @@ colvar::h_bond::h_bond(std::string const &conf) x.type(colvarvalue::type_scalar); init_scalar_boundaries(0.0, 1.0); - colvarproxy *proxy = cvm::main()->proxy; - int a_num = -1, d_num = -1; get_keyval(conf, "acceptor", a_num, a_num); get_keyval(conf, "donor", d_num, a_num); if ( (a_num == -1) || (d_num == -1) ) { - cvm::error("Error: either acceptor or donor undefined.\n"); - return; + error_code |= cvm::error("Error: either acceptor or donor undefined.\n", COLVARS_INPUT_ERROR); } cvm::atom acceptor = cvm::atom(a_num); @@ -327,34 +323,34 @@ colvar::h_bond::h_bond(std::string const &conf) atom_groups[0]->add_atom(acceptor); atom_groups[0]->add_atom(donor); - get_keyval(conf, "cutoff", r0, proxy->angstrom_to_internal(3.3)); - get_keyval(conf, "expNumer", en, 6); - get_keyval(conf, "expDenom", ed, 8); + get_keyval(conf, "cutoff", r0, r0); + get_keyval(conf, "expNumer", en, en); + get_keyval(conf, "expDenom", ed, ed); - if ( (en%2) || (ed%2) ) { - cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); + if ((en % 2) || (ed % 2)) { + error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", + COLVARS_INPUT_ERROR); } - if ( (en <= 0) || (ed <= 0) ) { - cvm::error("Error: negative exponent(s) provided.\n", - COLVARS_INPUT_ERROR); + if ((en <= 0) || (ed <= 0)) { + error_code |= cvm::error("Error: negative exponent(s) provided.\n", COLVARS_INPUT_ERROR); } if (cvm::debug()) cvm::log("Done initializing h_bond object.\n"); + + return error_code; } colvar::h_bond::h_bond(cvm::atom const &acceptor, cvm::atom const &donor, cvm::real r0_i, int en_i, int ed_i) - : r0(r0_i), en(en_i), ed(ed_i) + : h_bond() { - set_function_type("hBond"); - x.type(colvarvalue::type_scalar); - init_scalar_boundaries(0.0, 1.0); - + r0 = r0_i; + en = en_i; + ed = ed_i; register_atom_group(new cvm::atom_group); atom_groups[0]->add_atom(acceptor); atom_groups[0]->add_atom(donor); @@ -384,64 +380,59 @@ void colvar::h_bond::calc_gradients() } -void colvar::h_bond::apply_force(colvarvalue const &force) -{ - (atom_groups[0])->apply_colvar_force(force); -} - -simple_scalar_dist_functions(h_bond) - - - -colvar::selfcoordnum::selfcoordnum(std::string const &conf) - : cvc(conf), pairlist(NULL) +colvar::selfcoordnum::selfcoordnum() { set_function_type("selfCoordNum"); x.type(colvarvalue::type_scalar); + r0 = cvm::main()->proxy->angstrom_to_internal(4.0); +} - colvarproxy *proxy = cvm::main()->proxy; + +int colvar::selfcoordnum::init(std::string const &conf) +{ + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); - get_keyval(conf, "cutoff", r0, cvm::real(proxy->angstrom_to_internal(4.0))); - get_keyval(conf, "expNumer", en, 6); - get_keyval(conf, "expDenom", ed, 12); + get_keyval(conf, "cutoff", r0, r0); + get_keyval(conf, "expNumer", en, en); + get_keyval(conf, "expDenom", ed, ed); - if ( (en%2) || (ed%2) ) { - cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); + if ((en % 2) || (ed % 2)) { + error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", + COLVARS_INPUT_ERROR); } - if ( (en <= 0) || (ed <= 0) ) { - cvm::error("Error: negative exponent(s) provided.\n", - COLVARS_INPUT_ERROR); + if ((en <= 0) || (ed <= 0)) { + error_code |= cvm::error("Error: negative exponent(s) provided.\n", COLVARS_INPUT_ERROR); } if (!is_enabled(f_cvc_pbc_minimum_image)) { cvm::log("Warning: only minimum-image distances are used by this variable.\n"); } - get_keyval(conf, "tolerance", tolerance, 0.0); + get_keyval(conf, "tolerance", tolerance, tolerance); if (tolerance > 0) { - get_keyval(conf, "pairListFrequency", pairlist_freq, 100); + get_keyval(conf, "pairListFrequency", pairlist_freq, pairlist_freq); if ( ! (pairlist_freq > 0) ) { - cvm::error("Error: non-positive pairlistfrequency provided.\n", - COLVARS_INPUT_ERROR); - return; + error_code |= cvm::error("Error: non-positive pairlistfrequency provided.\n", + COLVARS_INPUT_ERROR); } pairlist = new bool[(group1->size()-1) * (group1->size()-1)]; } init_scalar_boundaries(0.0, static_cast((group1->size()-1) * (group1->size()-1))); + + return error_code; } colvar::selfcoordnum::~selfcoordnum() { - if (pairlist != NULL) { + if (pairlist) { delete [] pairlist; } } @@ -526,44 +517,36 @@ void colvar::selfcoordnum::calc_gradients() } -void colvar::selfcoordnum::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) { - group1->apply_colvar_force(force.real_value); - } -} - - -simple_scalar_dist_functions(selfcoordnum) - - -// groupcoordnum member functions -colvar::groupcoordnum::groupcoordnum(std::string const &conf) - : distance(conf), b_anisotropic(false) +colvar::groupcoordnum::groupcoordnum() { set_function_type("groupCoord"); x.type(colvarvalue::type_scalar); init_scalar_boundaries(0.0, 1.0); - colvarproxy *proxy = cvm::main()->proxy; + r0 = proxy->angstrom_to_internal(4.0); + r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0)); +} + + +int colvar::groupcoordnum::init(std::string const &conf) +{ + int error_code = distance::init(conf); // group1 and group2 are already initialized by distance() if (group1->b_dummy || group2->b_dummy) { - cvm::error("Error: neither group can be a dummy atom\n"); - return; + return cvm::error("Error: neither group can be a dummy atom\n", COLVARS_INPUT_ERROR); } - bool const b_scale = get_keyval(conf, "cutoff", r0, - cvm::real(proxy->angstrom_to_internal(4.0))); - - if (get_keyval(conf, "cutoff3", r0_vec, - cvm::rvector(4.0, 4.0, 4.0), parse_silent)) { + bool const b_scale = get_keyval(conf, "cutoff", r0, r0); + if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) { if (b_scale) { - cvm::error("Error: cannot specify \"scale\" and " - "\"scale3\" at the same time.\n"); - return; + error_code |= + cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n", + COLVARS_INPUT_ERROR); } b_anisotropic = true; // remove meaningless negative signs @@ -572,23 +555,23 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf) if (r0_vec.z < 0.0) r0_vec.z *= -1.0; } - get_keyval(conf, "expNumer", en, 6); - get_keyval(conf, "expDenom", ed, 12); + get_keyval(conf, "expNumer", en, en); + get_keyval(conf, "expDenom", ed, ed); - if ( (en%2) || (ed%2) ) { - cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); + if ((en % 2) || (ed % 2)) { + error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", + COLVARS_INPUT_ERROR); } - if ( (en <= 0) || (ed <= 0) ) { - cvm::error("Error: negative exponent(s) provided.\n", - COLVARS_INPUT_ERROR); + if ((en <= 0) || (ed <= 0)) { + error_code |= cvm::error("Error: negative exponent(s) provided.\n", COLVARS_INPUT_ERROR); } if (!is_enabled(f_cvc_pbc_minimum_image)) { cvm::log("Warning: only minimum-image distances are used by this variable.\n"); } + return error_code; } @@ -639,16 +622,3 @@ void colvar::groupcoordnum::calc_gradients() group1->set_weighted_gradient(group1_com_atom.grad); group2->set_weighted_gradient(group2_com_atom.grad); } - - -void colvar::groupcoordnum::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(groupcoordnum) diff --git a/src/colvarcomp_distances.cpp b/src/colvarcomp_distances.cpp index 589519e40..1ba6f99c3 100644 --- a/src/colvarcomp_distances.cpp +++ b/src/colvarcomp_distances.cpp @@ -16,8 +16,7 @@ #include "colvar_rotation_derivative.h" -colvar::distance::distance(std::string const &conf) - : cvc(conf) +colvar::distance::distance() { set_function_type("distance"); init_as_distance(); @@ -25,23 +24,19 @@ colvar::distance::distance(std::string const &conf) provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); - - group1 = parse_group(conf, "group1"); - group2 = parse_group(conf, "group2"); - - init_total_force_params(conf); } -colvar::distance::distance() - : cvc() +int colvar::distance::init(std::string const &conf) { - set_function_type("distance"); - init_as_distance(); + int error_code = cvc::init(conf); - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); + + error_code |= init_total_force_params(conf); + + return error_code; } @@ -83,35 +78,10 @@ void colvar::distance::calc_Jacobian_derivative() } -void colvar::distance::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(distance) - - - -colvar::distance_vec::distance_vec(std::string const &conf) - : distance(conf) -{ - set_function_type("distanceVec"); - enable(f_cvc_com_based); - disable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_3vector); -} - colvar::distance_vec::distance_vec() - : distance() { set_function_type("distanceVec"); - enable(f_cvc_com_based); disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_3vector); } @@ -145,48 +115,47 @@ void colvar::distance_vec::apply_force(colvarvalue const &force) } -cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, - colvarvalue const &x2) const +cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, colvarvalue const &x2) const { - return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2(); + if (is_enabled(f_cvc_pbc_minimum_image)) { + return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2(); + } + return (x2.rvector_value - x1.rvector_value).norm2(); } -colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); + if (is_enabled(f_cvc_pbc_minimum_image)) { + return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); + } + return 2.0 * (x2.rvector_value - x1.rvector_value); } -colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const +colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { - return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); + return distance_vec::dist2_lgrad(x2, x1); } +void colvar::distance_vec::wrap(colvarvalue & /* x_unwrapped */) const {} + -colvar::distance_z::distance_z(std::string const &conf) - : cvc(conf) +colvar::distance_z::distance_z() { set_function_type("distanceZ"); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); + provide(f_cvc_periodic); x.type(colvarvalue::type_scalar); +} - // TODO detect PBC from MD engine (in simple cases) - // and then update period in real time - if (period != 0.0) { - enable(f_cvc_periodic); - } - if ((wrap_center != 0.0) && !is_enabled(f_cvc_periodic)) { - cvm::error("Error: wrapAround was defined in a distanceZ component," - " but its period has not been set.\n"); - return; - } +int colvar::distance_z::init(std::string const &conf) +{ + int error_code = cvc::init(conf); main = parse_group(conf, "main"); ref1 = parse_group(conf, "ref"); @@ -201,8 +170,7 @@ colvar::distance_z::distance_z(std::string const &conf) } else { if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) { - cvm::error("Axis vector is zero!"); - return; + error_code |= cvm::error("Axis vector is zero!", COLVARS_INPUT_ERROR); } if (axis.norm2() != 1.0) { axis = axis.unit(); @@ -212,18 +180,9 @@ colvar::distance_z::distance_z(std::string const &conf) fixed_axis = true; } - init_total_force_params(conf); - -} - + error_code |= init_total_force_params(conf); -colvar::distance_z::distance_z() -{ - set_function_type("distanceZ"); - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); - x.type(colvarvalue::type_scalar); + return error_code; } @@ -253,7 +212,7 @@ void colvar::distance_z::calc_value() axis = axis.unit(); } x.real_value = axis * dist_v; - this->wrap(x); + wrap(x); } @@ -302,90 +261,12 @@ void colvar::distance_z::calc_Jacobian_derivative() } -void colvar::distance_z::apply_force(colvarvalue const &force) -{ - if (!ref1->noforce) - ref1->apply_colvar_force(force.real_value); - - if (ref2 && !ref2->noforce) - ref2->apply_colvar_force(force.real_value); - - if (!main->noforce) - main->apply_colvar_force(force.real_value); -} - - -// Differences should always be wrapped around 0 (ignoring wrap_center) -cvm::real colvar::distance_z::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - if (is_enabled(f_cvc_periodic)) { - cvm::real shift = cvm::floor(diff/period + 0.5); - diff -= shift * period; - } - return diff * diff; -} - - -colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - if (is_enabled(f_cvc_periodic)) { - cvm::real shift = cvm::floor(diff/period + 0.5); - diff -= shift * period; - } - return 2.0 * diff; -} - - -colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - if (is_enabled(f_cvc_periodic)) { - cvm::real shift = cvm::floor(diff/period + 0.5); - diff -= shift * period; - } - return (-2.0) * diff; -} - - -void colvar::distance_z::wrap(colvarvalue &x_unwrapped) const -{ - if (!is_enabled(f_cvc_periodic)) { - // don't wrap if the period has not been set - return; - } - cvm::real shift = - cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5); - x_unwrapped.real_value -= shift * period; -} - - - -colvar::distance_xy::distance_xy(std::string const &conf) - : distance_z(conf) -{ - set_function_type("distanceXY"); - init_as_distance(); - - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); -} - colvar::distance_xy::distance_xy() - : distance_z() { set_function_type("distanceXY"); + provide(f_cvc_periodic, false); // Disable inherited distance_z flag init_as_distance(); - - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - enable(f_cvc_com_based); } @@ -461,35 +342,8 @@ void colvar::distance_xy::calc_Jacobian_derivative() } -void colvar::distance_xy::apply_force(colvarvalue const &force) -{ - if (!ref1->noforce) - ref1->apply_colvar_force(force.real_value); - - if (ref2 && !ref2->noforce) - ref2->apply_colvar_force(force.real_value); - - if (!main->noforce) - main->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(distance_xy) - - - -colvar::distance_dir::distance_dir(std::string const &conf) - : distance(conf) -{ - set_function_type("distanceDir"); - enable(f_cvc_com_based); - disable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_unit3vector); -} - colvar::distance_dir::distance_dir() - : distance() { set_function_type("distanceDir"); enable(f_cvc_com_based); @@ -555,31 +409,39 @@ colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, } +void colvar::distance_dir::wrap(colvarvalue & /* x_unwrapped */) const {} + -colvar::distance_inv::distance_inv(std::string const &conf) - : cvc(conf) + +colvar::distance_inv::distance_inv() { set_function_type("distanceInv"); init_as_distance(); +} + + +int colvar::distance_inv::init(std::string const &conf) +{ + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); - get_keyval(conf, "exponent", exponent, 6); - if (exponent%2) { - cvm::error("Error: odd exponent provided, can only use even ones.\n"); - return; + get_keyval(conf, "exponent", exponent, exponent); + if (exponent % 2) { + error_code |= + cvm::error("Error: odd exponent provided, can only use even ones.\n", COLVARS_INPUT_ERROR); } if (exponent <= 0) { - cvm::error("Error: negative or zero exponent provided.\n"); - return; + error_code |= cvm::error("Error: negative or zero exponent provided.\n", COLVARS_INPUT_ERROR); } for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { if (ai1->id == ai2->id) { - cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n"); - return; + error_code |= cvm::error("Error: group1 and group2 have some atoms in common: this is not " + "allowed for distanceInv.\n", + COLVARS_INPUT_ERROR); } } } @@ -589,6 +451,8 @@ colvar::distance_inv::distance_inv(std::string const &conf) "for distanceInv, because its value and gradients are computed " "simultaneously.\n"); } + + return error_code; } @@ -641,39 +505,24 @@ void colvar::distance_inv::calc_gradients() } -void colvar::distance_inv::apply_force(colvarvalue const &force) -{ - if (!group1->noforce) - group1->apply_colvar_force(force.real_value); - if (!group2->noforce) - group2->apply_colvar_force(force.real_value); +colvar::distance_pairs::distance_pairs() +{ + set_function_type("distancePairs"); + disable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_vector); } -simple_scalar_dist_functions(distance_inv) - - - -colvar::distance_pairs::distance_pairs(std::string const &conf) - : cvc(conf) +int colvar::distance_pairs::init(std::string const &conf) { - set_function_type("distancePairs"); + int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); - - x.type(colvarvalue::type_vector); - disable(f_cvc_explicit_gradient); x.vector1d_value.resize(group1->size() * group2->size()); -} - -colvar::distance_pairs::distance_pairs() -{ - set_function_type("distancePairs"); - disable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_vector); + return error_code; } @@ -739,26 +588,28 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force) } +cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return (x1.vector1d_value - x2.vector1d_value).norm2(); +} + -colvar::dipole_magnitude::dipole_magnitude(std::string const &conf) - : cvc(conf) +colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - set_function_type("dipoleMagnitude"); - atoms = parse_group(conf, "atoms"); - init_total_force_params(conf); - x.type(colvarvalue::type_scalar); + return 2.0 * (x1.vector1d_value - x2.vector1d_value); } -colvar::dipole_magnitude::dipole_magnitude(cvm::atom const &a1) +colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { - set_function_type("dipoleMagnitude"); - atoms = new cvm::atom_group(std::vector(1, a1)); - register_atom_group(atoms); - x.type(colvarvalue::type_scalar); + return distance_pairs::dist2_lgrad(x1, x2); } +void colvar::distance_pairs::wrap(colvarvalue & /* x_unwrapped */) const {} + + + colvar::dipole_magnitude::dipole_magnitude() { set_function_type("dipoleMagnitude"); @@ -766,6 +617,15 @@ colvar::dipole_magnitude::dipole_magnitude() } +int colvar::dipole_magnitude::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + atoms = parse_group(conf, "atoms"); + if (!atoms) error_code |= COLVARS_INPUT_ERROR; + return error_code; +} + + void colvar::dipole_magnitude::calc_value() { cvm::atom_pos const atomsCom = atoms->center_of_mass(); @@ -786,26 +646,20 @@ void colvar::dipole_magnitude::calc_gradients() } -void colvar::dipole_magnitude::apply_force(colvarvalue const &force) + +colvar::gyration::gyration() { - if (!atoms->noforce) { - atoms->apply_colvar_force(force.real_value); - } + set_function_type("gyration"); + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + init_as_distance(); } -simple_scalar_dist_functions(dipole_magnitude) - - - -colvar::gyration::gyration(std::string const &conf) - : cvc(conf) +int colvar::gyration::init(std::string const &conf) { - set_function_type("gyration"); - init_as_distance(); + int error_code = cvc::init(conf); - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); atoms = parse_group(conf, "atoms"); if (atoms->b_user_defined_fit) { @@ -815,6 +669,8 @@ colvar::gyration::gyration(std::string const &conf) atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0)); } + + return error_code; } @@ -856,22 +712,10 @@ void colvar::gyration::calc_Jacobian_derivative() } -void colvar::gyration::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(gyration) - - -colvar::inertia::inertia(std::string const &conf) - : gyration(conf) +colvar::inertia::inertia() { set_function_type("inertia"); - init_as_distance(); } @@ -892,32 +736,26 @@ void colvar::inertia::calc_gradients() } -void colvar::inertia::apply_force(colvarvalue const &force) + +colvar::inertia_z::inertia_z() { - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); + set_function_type("inertiaZ"); } -simple_scalar_dist_functions(inertia_z) - - - -colvar::inertia_z::inertia_z(std::string const &conf) - : inertia(conf) +int colvar::inertia_z::init(std::string const &conf) { - set_function_type("inertiaZ"); - init_as_distance(); + int error_code = inertia::init(conf); if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) { - cvm::error("Axis vector is zero!", COLVARS_INPUT_ERROR); - return; + error_code |= cvm::error("Axis vector is zero!", COLVARS_INPUT_ERROR); } if (axis.norm2() != 1.0) { axis = axis.unit(); cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n"); } } + return error_code; } @@ -939,32 +777,23 @@ void colvar::inertia_z::calc_gradients() } -void colvar::inertia_z::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - -simple_scalar_dist_functions(inertia) - - - -colvar::rmsd::rmsd(std::string const &conf) - : cvc(conf) +colvar::rmsd::rmsd() { set_function_type("rmsd"); init_as_distance(); - provide(f_cvc_inv_gradient); +} + + +int colvar::rmsd::init(std::string const &conf) +{ + int error_code = cvc::init(conf); atoms = parse_group(conf, "atoms"); - if (cvm::get_error()) return; if (!atoms || atoms->size() == 0) { - cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD."); - return; + return error_code | COLVARS_INPUT_ERROR; } bool b_Jacobian_derivative = true; @@ -981,20 +810,20 @@ colvar::rmsd::rmsd(std::string const &conf) if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { cvm::log("Using reference positions from configuration file to calculate the variable.\n"); if (ref_pos.size() != atoms->size()) { - cvm::error("Error: the number of reference positions provided ("+ - cvm::to_str(ref_pos.size())+ - ") does not match the number of atoms of group \"atoms\" ("+ - cvm::to_str(atoms->size())+").\n"); - return; + error_code |= cvm::error("Error: the number of reference positions provided (" + + cvm::to_str(ref_pos.size()) + + ") does not match the number of atoms of group \"atoms\" (" + + cvm::to_str(atoms->size()) + ").\n", + COLVARS_INPUT_ERROR); } } else { // Only look for ref pos file if ref positions not already provided std::string ref_pos_file; if (get_keyval(conf, "refPositionsFile", ref_pos_file, std::string(""))) { if (ref_pos.size()) { - cvm::error("Error: cannot specify \"refPositionsFile\" and " - "\"refPositions\" at the same time.\n"); - return; + error_code |= cvm::error("Error: cannot specify \"refPositionsFile\" and " + "\"refPositions\" at the same time.\n", + COLVARS_INPUT_ERROR); } std::string ref_pos_col; @@ -1004,9 +833,9 @@ colvar::rmsd::rmsd(std::string const &conf) // if provided, use PDB column to select coordinates bool found = get_keyval(conf, "refPositionsColValue", ref_pos_col_value, 0.0); if (found && ref_pos_col_value==0.0) { - cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); - return; + error_code |= cvm::error("Error: refPositionsColValue, " + "if provided, must be non-zero.\n", + COLVARS_INPUT_ERROR); } } @@ -1015,15 +844,17 @@ colvar::rmsd::rmsd(std::string const &conf) cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms, ref_pos_col, ref_pos_col_value); } else { - cvm::error("Error: no reference positions for RMSD; use either refPositions of refPositionsFile."); - return; + error_code |= cvm::error( + "Error: no reference positions for RMSD; use either refPositions of refPositionsFile.", + COLVARS_INPUT_ERROR); } } if (ref_pos.size() != atoms->size()) { - cvm::error("Error: found " + cvm::to_str(ref_pos.size()) + - " reference positions for RMSD; expected " + cvm::to_str(atoms->size())); - return; + error_code |= + cvm::error("Error: found " + cvm::to_str(ref_pos.size()) + + " reference positions for RMSD; expected " + cvm::to_str(atoms->size()), + COLVARS_INPUT_ERROR); } if (atoms->b_user_defined_fit) { @@ -1045,6 +876,15 @@ colvar::rmsd::rmsd(std::string const &conf) } atoms->setup_rotation_derivative(); + error_code |= init_permutation(conf); + + return error_code; +} + + +int colvar::rmsd::init_permutation(std::string const &conf) +{ + int error_code = COLVARS_OK; std::string perm_conf; size_t pos = 0; // current position in config string n_permutations = 1; @@ -1059,21 +899,22 @@ colvar::rmsd::rmsd(std::string const &conf) std::vector const &ids = atoms->ids(); size_t const ia = std::find(ids.begin(), ids.end(), index-1) - ids.begin(); if (ia == atoms->size()) { - cvm::error("Error: atom id " + cvm::to_str(index) + - " is not a member of group \"atoms\"."); - return; + error_code |= cvm::error("Error: atom id " + cvm::to_str(index) + + " is not a member of group \"atoms\".", + COLVARS_INPUT_ERROR); } if (std::find(perm.begin(), perm.end(), ia) != perm.end()) { - cvm::error("Error: atom id " + cvm::to_str(index) + - " is mentioned more than once in atomPermutation list."); - return; + error_code |= cvm::error("Error: atom id " + cvm::to_str(index) + + " is mentioned more than once in atomPermutation list.", + COLVARS_INPUT_ERROR); } perm.push_back(ia); } if (perm.size() != atoms->size()) { - cvm::error("Error: symmetry permutation in input contains " + cvm::to_str(perm.size()) + - " indices, but group \"atoms\" contains " + cvm::to_str(atoms->size()) + " atoms."); - return; + error_code |= cvm::error( + "Error: symmetry permutation in input contains " + cvm::to_str(perm.size()) + + " indices, but group \"atoms\" contains " + cvm::to_str(atoms->size()) + " atoms.", + COLVARS_INPUT_ERROR); } cvm::log("atomPermutation = " + cvm::to_str(perm)); n_permutations++; @@ -1083,6 +924,8 @@ colvar::rmsd::rmsd(std::string const &conf) } } } + + return error_code; } @@ -1127,13 +970,6 @@ void colvar::rmsd::calc_gradients() } -void colvar::rmsd::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - void colvar::rmsd::calc_force_invgrads() { atoms->read_total_forces(); @@ -1210,20 +1046,24 @@ void colvar::rmsd::calc_Jacobian_derivative() } -simple_scalar_dist_functions(rmsd) - - -colvar::eigenvector::eigenvector(std::string const &conf) - : cvc(conf) +colvar::eigenvector::eigenvector() { set_function_type("eigenvector"); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); x.type(colvarvalue::type_scalar); +} + + +int colvar::eigenvector::init(std::string const &conf) +{ + int error_code = cvc::init(conf); atoms = parse_group(conf, "atoms"); - if (cvm::get_error()) return; + if (!atoms || atoms->size() == 0) { + return error_code | COLVARS_INPUT_ERROR; + } { bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos); @@ -1231,9 +1071,9 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (b_inline) { cvm::log("Using reference positions from input file.\n"); if (ref_pos.size() != atoms->size()) { - cvm::error("Error: reference positions do not " - "match the number of requested atoms.\n"); - return; + error_code |= cvm::error("Error: reference positions do not " + "match the number of requested atoms.\n", + COLVARS_INPUT_ERROR); } } @@ -1241,19 +1081,18 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (get_keyval(conf, "refPositionsFile", file_name)) { if (b_inline) { - cvm::error("Error: refPositions and refPositionsFile cannot be specified at the same time.\n"); - return; + error_code |= cvm::error( + "Error: refPositions and refPositionsFile cannot be specified at the same time.\n", + COLVARS_INPUT_ERROR); } - std::string file_col; double file_col_value=0.0; if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) { // use PDB flags if column is provided bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0); - if (found && file_col_value==0.0) { - cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); - return; + if (found && file_col_value == 0.0) { + error_code |= cvm::error("Error: refPositionsColValue, if provided, must be non-zero.\n", + COLVARS_INPUT_ERROR); } } @@ -1264,14 +1103,14 @@ colvar::eigenvector::eigenvector(std::string const &conf) } if (ref_pos.size() == 0) { - cvm::error("Error: reference positions were not provided.\n", COLVARS_INPUT_ERROR); - return; + error_code |= + cvm::error("Error: reference positions were not provided.\n", COLVARS_INPUT_ERROR); } if (ref_pos.size() != atoms->size()) { - cvm::error("Error: reference positions do not " - "match the number of requested atoms.\n", COLVARS_INPUT_ERROR); - return; + error_code |= cvm::error("Error: reference positions do not " + "match the number of requested atoms.\n", + COLVARS_INPUT_ERROR); } // save for later the geometric center of the provided positions (may not be the origin) @@ -1302,9 +1141,8 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (b_inline) { cvm::log("Using vector components from input file.\n"); if (eigenvec.size() != atoms->size()) { - cvm::error("Error: vector components do not " - "match the number of requested atoms->\n"); - return; + error_code |= cvm::error("Error: vector components do not " + "match the number of requested atoms->\n", COLVARS_INPUT_ERROR); } } @@ -1312,8 +1150,9 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (get_keyval(conf, "vectorFile", file_name)) { if (b_inline) { - cvm::error("Error: vector and vectorFile cannot be specified at the same time.\n"); - return; + error_code |= + cvm::error("Error: vector and vectorFile cannot be specified at the same time.\n", + COLVARS_INPUT_ERROR); } std::string file_col; @@ -1322,8 +1161,8 @@ colvar::eigenvector::eigenvector(std::string const &conf) // use PDB flags if column is provided bool found = get_keyval(conf, "vectorColValue", file_col_value, 0.0); if (found && file_col_value==0.0) { - cvm::error("Error: vectorColValue, if provided, must be non-zero.\n"); - return; + error_code |= cvm::error("Error: vectorColValue, if provided, must be non-zero.\n", + COLVARS_INPUT_ERROR); } } @@ -1334,9 +1173,8 @@ colvar::eigenvector::eigenvector(std::string const &conf) } if (!ref_pos.size() || !eigenvec.size()) { - cvm::error("Error: both reference coordinates" - "and eigenvector must be defined.\n"); - return; + error_code |= cvm::error("Error: both reference coordinates and eigenvector must be defined.\n", + COLVARS_INPUT_ERROR); } cvm::atom_pos eig_center(0.0, 0.0, 0.0); @@ -1410,6 +1248,8 @@ colvar::eigenvector::eigenvector(std::string const &conf) cvm::log("The norm of the vector is |v| = "+ cvm::to_str(1.0/cvm::sqrt(eigenvec_invnorm2))+".\n"); } + + return error_code; } @@ -1430,13 +1270,6 @@ void colvar::eigenvector::calc_gradients() } -void colvar::eigenvector::apply_force(colvarvalue const &force) -{ - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); -} - - void colvar::eigenvector::calc_force_invgrads() { atoms->read_total_forces(); @@ -1502,14 +1335,18 @@ void colvar::eigenvector::calc_Jacobian_derivative() } -simple_scalar_dist_functions(eigenvector) +colvar::cartesian::cartesian() +{ + set_function_type("cartesian"); + x.type(colvarvalue::type_vector); + disable(f_cvc_explicit_gradient); +} -colvar::cartesian::cartesian(std::string const &conf) - : cvc(conf) +int colvar::cartesian::init(std::string const &conf) { - set_function_type("cartesian"); + int error_code = cvc::init(conf); atoms = parse_group(conf, "atoms"); @@ -1524,14 +1361,15 @@ colvar::cartesian::cartesian(std::string const &conf) if (use_z) axes.push_back(2); if (axes.size() == 0) { - cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n"); - return; + error_code |= + cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n", + COLVARS_INPUT_ERROR); } - x.type(colvarvalue::type_vector); - disable(f_cvc_explicit_gradient); // Don't try to access atoms if creation of the atom group failed - if (atoms != NULL) x.vector1d_value.resize(atoms->size() * axes.size()); + if (atoms) x.vector1d_value.resize(atoms->size() * axes.size()); + + return error_code; } @@ -1569,3 +1407,24 @@ void colvar::cartesian::apply_force(colvarvalue const &force) } } } + + +cvm::real colvar::cartesian::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return (x1.vector1d_value - x2.vector1d_value).norm2(); +} + + +colvarvalue colvar::cartesian::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return 2.0 * (x1.vector1d_value - x2.vector1d_value); +} + + +colvarvalue colvar::cartesian::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return cartesian::dist2_lgrad(x1, x2); +} + + +void colvar::cartesian::wrap(colvarvalue & /* x_unwrapped */) const {} diff --git a/src/colvarcomp_gpath.cpp b/src/colvarcomp_gpath.cpp index 0b9a18c43..46d6b8e5a 100644 --- a/src/colvarcomp_gpath.cpp +++ b/src/colvarcomp_gpath.cpp @@ -18,7 +18,21 @@ #include "colvar.h" #include "colvarcomp.h" -colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(conf), atoms(nullptr), reference_frames(0) { + + +colvar::CartesianBasedPath::CartesianBasedPath() +{ + x.type(colvarvalue::type_scalar); + // Don't use implicit gradient + enable(f_cvc_explicit_gradient); +} + + +int colvar::CartesianBasedPath::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + if (error_code != COLVARS_OK) return error_code; + // Parse selected atoms atoms = parse_group(conf, "atoms"); has_user_defined_fitting = false; @@ -33,9 +47,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con if (get_keyval(conf, "refPositionsCol", reference_column, std::string(""))) { bool found = get_keyval(conf, "refPositionsColValue", reference_column_value, reference_column_value); if (found && reference_column_value == 0.0) { - cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); - return; + return cvm::error("Error: refPositionsColValue, if provided, must be non-zero.\n", + COLVARS_INPUT_ERROR); } } // Lookup all reference frames @@ -91,9 +104,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con tmp_atoms->setup_rotation_derivative(); comp_atoms.push_back(tmp_atoms); } - x.type(colvarvalue::type_scalar); - // Don't use implicit gradient - enable(f_cvc_explicit_gradient); + + return error_code; } colvar::CartesianBasedPath::~CartesianBasedPath() { @@ -153,8 +165,26 @@ void colvar::CartesianBasedPath::computeDistanceBetweenReferenceFrames(std::vect } } -colvar::gspath::gspath(std::string const &conf): CartesianBasedPath(conf) { + +void colvar::CartesianBasedPath::apply_force(colvarvalue const &force) +{ + cvm::error("Error: using apply_force() in a component of type CartesianBasedPath, which is abstract.\n", + COLVARS_BUG_ERROR); +} + + + +colvar::gspath::gspath() +{ set_function_type("gspath"); +} + + +int colvar::gspath::init(std::string const &conf) +{ + int error_code = CartesianBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; + get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true); if (use_second_closest_frame == true) { cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n")); @@ -168,12 +198,14 @@ colvar::gspath::gspath(std::string const &conf): CartesianBasedPath(conf) { cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n")); } if (total_reference_frames < 2) { - cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspath requires at least 2 frames.\n"); - return; + return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + + " reference frames, but gspath requires at least 2 frames.\n", + COLVARS_INPUT_ERROR); } GeometricPathCV::GeometricPathBase::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame); cvm::log(std::string("Geometric pathCV(s) is initialized.\n")); cvm::log(std::string("Geometric pathCV(s) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n")); + return error_code; } void colvar::gspath::updateDistanceToReferenceFrames() { @@ -297,8 +329,17 @@ void colvar::gspath::apply_force(colvarvalue const &force) { (*(comp_atoms[min_frame_index_2])).apply_colvar_force(F); } -colvar::gzpath::gzpath(std::string const &conf): CartesianBasedPath(conf) { + +colvar::gzpath::gzpath() +{ set_function_type("gzpath"); +} + +int colvar::gzpath::init(std::string const &conf) +{ + int error_code = CartesianBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; + get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true); if (use_second_closest_frame == true) { cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n")); @@ -317,13 +358,15 @@ colvar::gzpath::gzpath(std::string const &conf): CartesianBasedPath(conf) { cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n")); } if (total_reference_frames < 2) { - cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpath requires at least 2 frames.\n"); - return; + return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + + " reference frames, but gzpath requires at least 2 frames.\n", + COLVARS_INPUT_ERROR); } GeometricPathCV::GeometricPathBase::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square); // Logging cvm::log(std::string("Geometric pathCV(z) is initialized.\n")); cvm::log(std::string("Geometric pathCV(z) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n")); + return error_code; } void colvar::gzpath::updateDistanceToReferenceFrames() { @@ -431,14 +474,26 @@ void colvar::gzpath::apply_force(colvarvalue const &force) { } -colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { +colvar::CVBasedPath::CVBasedPath() +{ + set_function_type("gspathCV"); + x.type(colvarvalue::type_scalar); +} + + +int colvar::CVBasedPath::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + if (error_code != COLVARS_OK) return error_code; + // Lookup all available sub-cvcs for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) { if (key_lookup(conf, it_cv_map->first.c_str())) { std::vector sub_cvc_confs; get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs); for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) { - cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf))); + cv.push_back((it_cv_map->second)()); + cv.back()->init(*(it_sub_cvc_conf)); } } } @@ -461,7 +516,7 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { cvm::log(std::string("Reading path file: ") + path_filename + std::string("\n")); auto &ifs_path = cvm::main()->proxy->input_stream(path_filename); if (!ifs_path) { - return; + return COLVARS_INPUT_ERROR; } std::string line; const std::string token(" "); @@ -482,8 +537,8 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { cvm::log(cvm::to_str(tmp_cv[i_cv][i - start_index])); } } else { - cvm::error("Error: incorrect format of path file.\n"); - return; + error_code = cvm::error("Error: incorrect format of path file.\n", COLVARS_INPUT_ERROR); + return error_code; } } if (!fields.empty()) { @@ -493,15 +548,18 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { } cvm::main()->proxy->close_input_stream(path_filename); if (total_reference_frames <= 1) { - cvm::error("Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n"); - return; + error_code = cvm::error( + "Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n", + COLVARS_INPUT_ERROR); + return error_code; } if (cv.size() == 0) { - cvm::error("Error: the CV " + name + - " expects one or more nesting components.\n"); - return; + error_code = + cvm::error("Error: the CV " + name + " expects one or more nesting components.\n", + COLVARS_INPUT_ERROR); + return error_code; } - x.type(colvarvalue::type_scalar); + use_explicit_gradients = true; for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { @@ -511,6 +569,8 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { if (!use_explicit_gradients) { disable(f_cvc_explicit_gradient); } + + return error_code; } void colvar::CVBasedPath::computeDistanceToReferenceFrames(std::vector& result) { @@ -580,8 +640,47 @@ colvar::CVBasedPath::~CVBasedPath() { atom_groups.clear(); } -colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) { + +void colvar::CVBasedPath::apply_force(colvarvalue const &force) +{ + cvm::error("Error: using apply_force() in a component of type CVBasedPath, which is abstract.\n", + COLVARS_BUG_ERROR); +} + + +cvm::real colvar::CVBasedPath::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x1.dist2(x2); +} + + +colvarvalue colvar::CVBasedPath::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x1.dist2_grad(x2); +} + + +colvarvalue colvar::CVBasedPath::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x2.dist2_grad(x1); +} + + +void colvar::CVBasedPath::wrap(colvarvalue & /* x_unwrapped */) const {} + + + +colvar::gspathCV::gspathCV() +{ set_function_type("gspathCV"); + x.type(colvarvalue::type_scalar); +} + +int colvar::gspathCV::init(std::string const &conf) +{ + int error_code = CVBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; + cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); // Initialize variables for future calculation get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true); @@ -597,11 +696,10 @@ colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) { cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n")); } if (total_reference_frames < 2) { - cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspathCV requires at least 2 frames.\n"); - return; + return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspathCV requires at least 2 frames.\n", COLVARS_INPUT_ERROR); } GeometricPathCV::GeometricPathBase::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame); - x.type(colvarvalue::type_scalar); + return error_code; } colvar::gspathCV::~gspathCV() {} @@ -703,8 +801,16 @@ void colvar::gspathCV::apply_force(colvarvalue const &force) { } } -colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) { +colvar::gzpathCV::gzpathCV() +{ set_function_type("gzpathCV"); +} + +int colvar::gzpathCV::init(std::string const &conf) { + + int error_code = CVBasedPath::init(conf); + if (error_code != COLVARS_OK) return error_code; + cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); // Initialize variables for future calculation M = cvm::real(total_reference_frames - 1); @@ -727,11 +833,13 @@ colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) { cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n")); } if (total_reference_frames < 2) { - cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpathCV requires at least 2 frames.\n"); - return; + return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + + " reference frames, but gzpathCV requires at least 2 frames.\n", + COLVARS_INPUT_ERROR); } GeometricPathCV::GeometricPathBase::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square); - x.type(colvarvalue::type_scalar); + + return error_code; } colvar::gzpathCV::~gzpathCV() { diff --git a/src/colvarcomp_neuralnetwork.cpp b/src/colvarcomp_neuralnetwork.cpp index e3550a5af..2d2f89d5f 100644 --- a/src/colvarcomp_neuralnetwork.cpp +++ b/src/colvarcomp_neuralnetwork.cpp @@ -15,8 +15,17 @@ using namespace neuralnetworkCV; -colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination(conf) { + +colvar::neuralNetwork::neuralNetwork() +{ set_function_type("neuralNetwork"); +} + + +int colvar::neuralNetwork::init(std::string const &conf) +{ + int error_code = linearCombination::init(conf); + if (error_code != COLVARS_OK) return error_code; // the output of neural network consists of multiple values // read "output_component" key to determine it get_keyval(conf, "output_component", m_output_index); @@ -65,8 +74,8 @@ colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination std::string function_name; get_keyval(conf, lookup_key.c_str(), function_name, std::string("")); if (activation_function_map.find(function_name) == activation_function_map.end()) { - cvm::error("Unknown activation function name: \"" + function_name + "\".\n"); - return; + return cvm::error("Unknown activation function name: \"" + function_name + "\".\n", + COLVARS_INPUT_ERROR); } activation_functions.push_back(std::make_pair(false, function_name)); cvm::log(std::string{"The activation function for layer["} + cvm::to_str(num_activation_functions + 1) + std::string{"] is "} + function_name + '\n'); @@ -85,11 +94,13 @@ colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination } // expect the three numbers are equal if ((num_layers_weight != num_layers_bias) || (num_layers_bias != num_activation_functions)) { - cvm::error("Error: the numbers of weights, biases and activation functions do not match.\n"); - return; + return cvm::error( + "Error: the numbers of weights, biases and activation functions do not match.\n", + COLVARS_INPUT_ERROR); } // nn = std::make_unique(); // std::make_unique is only available in C++14 + if (nn) nn.reset(); nn = std::unique_ptr(new neuralnetworkCV::neuralNetworkCompute()); for (size_t i_layer = 0; i_layer < num_layers_weight; ++i_layer) { denseLayer d; @@ -99,8 +110,9 @@ colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination try { d = denseLayer(weight_files[i_layer], bias_files[i_layer], activation_functions[i_layer].second); } catch (std::exception &ex) { - cvm::error("Error on initializing layer " + cvm::to_str(i_layer) + " (" + ex.what() + ")\n", COLVARS_INPUT_ERROR); - return; + return cvm::error("Error on initializing layer " + cvm::to_str(i_layer) + + " (" + ex.what() + ")\n", + COLVARS_INPUT_ERROR); } } else { #endif @@ -110,8 +122,9 @@ colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination try { d = denseLayer(weight_files[i_layer], bias_files[i_layer], f, df); } catch (std::exception &ex) { - cvm::error("Error on initializing layer " + cvm::to_str(i_layer) + " (" + ex.what() + ")\n", COLVARS_INPUT_ERROR); - return; + return cvm::error("Error on initializing layer " + cvm::to_str(i_layer) + + " (" + ex.what() + ")\n", + COLVARS_INPUT_ERROR); } #ifdef LEPTON } @@ -129,11 +142,11 @@ colvar::neuralNetwork::neuralNetwork(std::string const &conf): linearCombination } } } else { - cvm::error("Error: error on adding a new dense layer.\n"); - return; + return cvm::error("Error: error on adding a new dense layer.\n", COLVARS_INPUT_ERROR); } } nn->input().resize(cv.size()); + return error_code; } colvar::neuralNetwork::~neuralNetwork() { @@ -190,3 +203,25 @@ void colvar::neuralNetwork::apply_force(colvarvalue const &force) { } } } + + +cvm::real colvar::neuralNetwork::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x1.dist2(x2); +} + + +colvarvalue colvar::neuralNetwork::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x1.dist2_grad(x2); +} + + +colvarvalue colvar::neuralNetwork::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const +{ + return x2.dist2_grad(x1); +} + + + +void colvar::neuralNetwork::wrap(colvarvalue & /* x_unwrapped */) const {} diff --git a/src/colvarcomp_protein.cpp b/src/colvarcomp_protein.cpp index 0f3830614..f78209514 100644 --- a/src/colvarcomp_protein.cpp +++ b/src/colvarcomp_protein.cpp @@ -15,14 +15,19 @@ #include "colvarcomp.h" -colvar::alpha_angles::alpha_angles(std::string const &conf) - : cvc(conf) +colvar::alpha_angles::alpha_angles() { set_function_type("alpha"); enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); - colvarproxy *proxy = cvm::main()->proxy; + r0 = proxy->angstrom_to_internal(3.3); +} + + +int colvar::alpha_angles::init(std::string const &conf) +{ + int error_code = cvc::init(conf); std::string segment_id; get_keyval(conf, "psfSegID", segment_id, std::string("MAIN")); @@ -43,29 +48,29 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) } } } else { - cvm::error("Error: no residues defined in \"residueRange\".\n"); - return; + error_code |= + cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR); } } if (residues.size() < 5) { - cvm::error("Error: not enough residues defined in \"residueRange\".\n"); - return; + error_code |= cvm::error("Error: not enough residues defined in \"residueRange\".\n", + COLVARS_INPUT_ERROR); } std::string const &sid = segment_id; std::vector const &r = residues; - get_keyval(conf, "hBondCoeff", hb_coeff, 0.5); - if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) { - cvm::error("Error: hBondCoeff must be defined between 0 and 1.\n"); - return; + get_keyval(conf, "hBondCoeff", hb_coeff, hb_coeff); + if ((hb_coeff < 0.0) || (hb_coeff > 1.0)) { + error_code |= + cvm::error("Error: hBondCoeff must be defined between 0 and 1.\n", COLVARS_INPUT_ERROR); } - get_keyval(conf, "angleRef", theta_ref, 88.0); - get_keyval(conf, "angleTol", theta_tol, 15.0); + get_keyval(conf, "angleRef", theta_ref, theta_ref); + get_keyval(conf, "angleTol", theta_tol, theta_tol); if (hb_coeff < 1.0) { @@ -83,11 +88,9 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) } { - cvm::real r0; - size_t en, ed; - get_keyval(conf, "hBondCutoff", r0, proxy->angstrom_to_internal(3.3)); - get_keyval(conf, "hBondExpNumer", en, 6); - get_keyval(conf, "hBondExpDenom", ed, 8); + get_keyval(conf, "hBondCutoff", r0, r0); + get_keyval(conf, "hBondExpNumer", en, en); + get_keyval(conf, "hBondExpDenom", ed, ed); if (hb_coeff > 0.0) { @@ -102,15 +105,8 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) cvm::log("The hBondCoeff specified will disable the hydrogen bond terms.\n"); } } -} - -colvar::alpha_angles::alpha_angles() - : cvc() -{ - set_function_type("alphaAngles"); - enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); + return error_code; } @@ -272,24 +268,27 @@ void colvar::alpha_angles::apply_force(colvarvalue const &force) } -simple_scalar_dist_functions(alpha_angles) - ////////////////////////////////////////////////////////////////////// // dihedral principal component ////////////////////////////////////////////////////////////////////// -colvar::dihedPC::dihedPC(std::string const &conf) - : cvc(conf) +colvar::dihedPC::dihedPC() { - if (cvm::debug()) - cvm::log("Initializing dihedral PC object.\n"); - set_function_type("dihedPC"); // Supported through references to atom groups of children cvcs enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); +} + + +int colvar::dihedPC::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + + if (cvm::debug()) + cvm::log("Initializing dihedral PC object.\n"); std::string segment_id; get_keyval(conf, "psfSegID", segment_id, std::string("MAIN")); @@ -310,14 +309,14 @@ colvar::dihedPC::dihedPC(std::string const &conf) } } } else { - cvm::error("Error: no residues defined in \"residueRange\".\n"); - return; + error_code |= + cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR); } } if (residues.size() < 2) { - cvm::error("Error: dihedralPC requires at least two residues.\n"); - return; + error_code |= + cvm::error("Error: dihedralPC requires at least two residues.\n", COLVARS_INPUT_ERROR); } std::string const &sid = segment_id; @@ -328,15 +327,15 @@ colvar::dihedPC::dihedPC(std::string const &conf) if (get_keyval(conf, "vectorFile", vecFileName, vecFileName)) { get_keyval(conf, "vectorNumber", vecNumber, 0); if (vecNumber < 1) { - cvm::error("A positive value of vectorNumber is required."); - return; + error_code |= + cvm::error("A positive value of vectorNumber is required.", COLVARS_INPUT_ERROR); } std::istream &vecFile = cvm::main()->proxy->input_stream(vecFileName, "dihedral PCA vector file"); if (!vecFile) { - return; + return COLVARS_INPUT_ERROR; } // TODO: adapt to different formats by setting this flag @@ -382,11 +381,10 @@ colvar::dihedPC::dihedPC(std::string const &conf) } if ( coeffs.size() != 4 * (residues.size() - 1)) { - cvm::error("Error: wrong number of coefficients: " + - cvm::to_str(coeffs.size()) + ". Expected " + - cvm::to_str(4 * (residues.size() - 1)) + - " (4 coeffs per residue, minus one residue).\n"); - return; + error_code |= cvm::error("Error: wrong number of coefficients: " + cvm::to_str(coeffs.size()) + + ". Expected " + cvm::to_str(4 * (residues.size() - 1)) + + " (4 coeffs per residue, minus one residue).\n", + COLVARS_INPUT_ERROR); } for (size_t i = 0; i < residues.size()-1; i++) { @@ -412,16 +410,8 @@ colvar::dihedPC::dihedPC(std::string const &conf) if (cvm::debug()) cvm::log("Done initializing dihedPC object.\n"); -} - -colvar::dihedPC::dihedPC() - : cvc() -{ - set_function_type("dihedPC"); - // Supported through references to atom groups of children cvcs - enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); + return error_code; } @@ -490,6 +480,3 @@ void colvar::dihedPC::apply_force(colvarvalue const &force) coeffs[2*i+1] * dsindt) * force); } } - - -simple_scalar_dist_functions(dihedPC) diff --git a/src/colvarcomp_rotations.cpp b/src/colvarcomp_rotations.cpp index d59f2c93e..d4aa9c024 100644 --- a/src/colvarcomp_rotations.cpp +++ b/src/colvarcomp_rotations.cpp @@ -22,13 +22,12 @@ struct colvar::orientation::rotation_derivative_impl_: public rotation_derivativ }; -colvar::orientation::orientation(std::string const &conf) - : cvc() +colvar::orientation::orientation() { set_function_type("orientation"); + rot_deriv_impl = std::unique_ptr(new rotation_derivative_impl_(this)); disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_quaternion); - colvar::orientation::init(conf); } @@ -94,8 +93,6 @@ int colvar::orientation::init(std::string const &conf) get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0)); - rot_deriv_impl = std::unique_ptr(new rotation_derivative_impl_(this)); - // If the debug gradients feature is active, debug the rotation gradients // (note that this won't be active for the orientation CVC itself, because // colvardeps prevents the flag's activation) @@ -105,15 +102,6 @@ int colvar::orientation::init(std::string const &conf) } -colvar::orientation::orientation() - : cvc() -{ - set_function_type("orientation"); - disable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_quaternion); -} - - void colvar::orientation::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -176,20 +164,15 @@ colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1, } +void colvar::orientation::wrap(colvarvalue & /* x_unwrapped */) const {} + + -colvar::orientation_angle::orientation_angle(std::string const &conf) - : orientation() +colvar::orientation_angle::orientation_angle() { set_function_type("orientationAngle"); init_as_angle(); enable(f_cvc_explicit_gradient); - orientation_angle::init(conf); -} - - -int colvar::orientation_angle::init(std::string const &conf) -{ - return orientation::init(conf); } @@ -226,31 +209,40 @@ void colvar::orientation_angle::calc_gradients() void colvar::orientation_angle::apply_force(colvarvalue const &force) { - cvm::real const &fw = force.real_value; - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } + cvc::apply_force(force); +} + + +cvm::real colvar::orientation_angle::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + return cvc::dist2(x1, x2); +} + + +colvarvalue colvar::orientation_angle::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return cvc::dist2_lgrad(x1, x2); +} + + +colvarvalue colvar::orientation_angle::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return cvc::dist2_rgrad(x1, x2); } -simple_scalar_dist_functions(orientation_angle) +void colvar::orientation_angle::wrap(colvarvalue & /* x_unwrapped */) const {} -colvar::orientation_proj::orientation_proj(std::string const &conf) - : orientation() +colvar::orientation_proj::orientation_proj() { set_function_type("orientationProj"); enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); init_scalar_boundaries(0.0, 1.0); - orientation_proj::init(conf); -} - - -int colvar::orientation_proj::init(std::string const &conf) -{ - return orientation::init(conf); } @@ -275,36 +267,19 @@ void colvar::orientation_proj::calc_gradients() } -void colvar::orientation_proj::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -simple_scalar_dist_functions(orientation_proj) - - - -colvar::tilt::tilt(std::string const &conf) - : orientation() +colvar::tilt::tilt() { set_function_type("tilt"); x.type(colvarvalue::type_scalar); enable(f_cvc_explicit_gradient); init_scalar_boundaries(-1.0, 1.0); - tilt::init(conf); } int colvar::tilt::init(std::string const &conf) { - int error_code = COLVARS_OK; - - error_code |= orientation::init(conf); + int error_code = orientation_proj::init(conf); get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0)); if (axis.norm2() != 1.0) { @@ -343,55 +318,12 @@ void colvar::tilt::calc_gradients() } -void colvar::tilt::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -simple_scalar_dist_functions(tilt) - - - -colvar::spin_angle::spin_angle(std::string const &conf) - : orientation() -{ - set_function_type("spinAngle"); - init_as_periodic_angle(); - enable(f_cvc_periodic); - enable(f_cvc_explicit_gradient); - spin_angle::init(conf); -} - - -int colvar::spin_angle::init(std::string const &conf) -{ - int error_code = COLVARS_OK; - - error_code |= orientation::init(conf); - - get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0)); - if (axis.norm2() != 1.0) { - axis /= axis.norm(); - cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n"); - } - - return error_code; -} - colvar::spin_angle::spin_angle() - : orientation() { set_function_type("spinAngle"); - period = 360.0; - enable(f_cvc_periodic); + init_as_periodic_angle(); enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); } @@ -403,7 +335,7 @@ void colvar::spin_angle::calc_value() rot.calc_optimal_rotation(ref_pos, shifted_pos); x.real_value = rot.spin_angle(axis); - this->wrap(x); + wrap(x); } @@ -423,71 +355,8 @@ void colvar::spin_angle::calc_gradients() } -void colvar::spin_angle::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -cvm::real colvar::spin_angle::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return diff * diff; -} - - -colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return 2.0 * diff; -} - - -colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return (-2.0) * diff; -} - - -void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const -{ - if ((x_unwrapped.real_value - wrap_center) >= 180.0) { - x_unwrapped.real_value -= 360.0; - return; - } - - if ((x_unwrapped.real_value - wrap_center) < -180.0) { - x_unwrapped.real_value += 360.0; - return; - } - - return; -} - - -colvar::euler_phi::euler_phi(std::string const &conf) - : orientation() -{ - set_function_type("eulerPhi"); - init_as_periodic_angle(); - enable(f_cvc_explicit_gradient); - euler_phi::init(conf); -} - colvar::euler_phi::euler_phi() - : orientation() { set_function_type("eulerPhi"); init_as_periodic_angle(); @@ -495,14 +364,6 @@ colvar::euler_phi::euler_phi() } -int colvar::euler_phi::init(std::string const &conf) -{ - int error_code = COLVARS_OK; - error_code |= orientation::init(conf); - return error_code; -} - - void colvar::euler_phi::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -543,70 +404,8 @@ void colvar::euler_phi::calc_gradients() } -void colvar::euler_phi::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -cvm::real colvar::euler_phi::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return diff * diff; -} - - -colvarvalue colvar::euler_phi::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return 2.0 * diff; -} - - -colvarvalue colvar::euler_phi::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return (-2.0) * diff; -} - - -void colvar::euler_phi::wrap(colvarvalue &x_unwrapped) const -{ - if ((x_unwrapped.real_value - wrap_center) >= 180.0) { - x_unwrapped.real_value -= 360.0; - return; - } - - if ((x_unwrapped.real_value - wrap_center) < -180.0) { - x_unwrapped.real_value += 360.0; - return; - } - - return; -} - - -colvar::euler_psi::euler_psi(std::string const &conf) - : orientation() -{ - set_function_type("eulerPsi"); - init_as_periodic_angle(); - enable(f_cvc_explicit_gradient); - euler_psi::init(conf); -} - colvar::euler_psi::euler_psi() - : orientation() { set_function_type("eulerPsi"); init_as_periodic_angle(); @@ -614,14 +413,6 @@ colvar::euler_psi::euler_psi() } -int colvar::euler_psi::init(std::string const &conf) -{ - int error_code = COLVARS_OK; - error_code |= orientation::init(conf); - return error_code; -} - - void colvar::euler_psi::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -662,70 +453,8 @@ void colvar::euler_psi::calc_gradients() } -void colvar::euler_psi::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -cvm::real colvar::euler_psi::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return diff * diff; -} - - -colvarvalue colvar::euler_psi::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return 2.0 * diff; -} - - -colvarvalue colvar::euler_psi::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - cvm::real diff = x1.real_value - x2.real_value; - diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); - return (-2.0) * diff; -} - - -void colvar::euler_psi::wrap(colvarvalue &x_unwrapped) const -{ - if ((x_unwrapped.real_value - wrap_center) >= 180.0) { - x_unwrapped.real_value -= 360.0; - return; - } - - if ((x_unwrapped.real_value - wrap_center) < -180.0) { - x_unwrapped.real_value += 360.0; - return; - } - - return; -} - - -colvar::euler_theta::euler_theta(std::string const &conf) - : orientation() -{ - set_function_type("eulerTheta"); - init_as_angle(); - enable(f_cvc_explicit_gradient); - euler_theta::init(conf); -} - colvar::euler_theta::euler_theta() - : orientation() { set_function_type("eulerTheta"); init_as_angle(); @@ -733,14 +462,6 @@ colvar::euler_theta::euler_theta() } -int colvar::euler_theta::init(std::string const &conf) -{ - int error_code = COLVARS_OK; - error_code |= orientation::init(conf); - return error_code; -} - - void colvar::euler_theta::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -777,36 +498,3 @@ void colvar::euler_theta::calc_gradients() (dxdq3 * dq0_2[3]); } } - - -void colvar::euler_theta::apply_force(colvarvalue const &force) -{ - cvm::real const &fw = force.real_value; - if (!atoms->noforce) { - atoms->apply_colvar_force(fw); - } -} - - -cvm::real colvar::euler_theta::dist2(colvarvalue const &x1, - colvarvalue const &x2) const -{ - // theta angle is not periodic - return cvc::dist2(x1, x2); -} - - -colvarvalue colvar::euler_theta::dist2_lgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - // theta angle is not periodic - return cvc::dist2_lgrad(x1, x2); -} - - -colvarvalue colvar::euler_theta::dist2_rgrad(colvarvalue const &x1, - colvarvalue const &x2) const -{ - // theta angle is not periodic - return cvc::dist2_rgrad(x1, x2); -} diff --git a/src/colvarcomp_volmaps.cpp b/src/colvarcomp_volmaps.cpp index a95f15bcb..a15501ab6 100644 --- a/src/colvarcomp_volmaps.cpp +++ b/src/colvarcomp_volmaps.cpp @@ -15,28 +15,12 @@ colvar::map_total::map_total() - : cvc() { set_function_type("mapTotal"); - volmap_id = -1; - volmap_index = -1; - atoms = NULL; x.type(colvarvalue::type_scalar); } -colvar::map_total::map_total(std::string const &conf) - : cvc() // init() will take care of this -{ - set_function_type("mapTotal"); - volmap_id = -1; - volmap_index = -1; - atoms = NULL; - x.type(colvarvalue::type_scalar); - map_total::init(conf); -} - - int colvar::map_total::init(std::string const &conf) { int error_code = cvc::init(conf); @@ -49,12 +33,12 @@ int colvar::map_total::init(std::string const &conf) if ((volmap_name.size() > 0) && (volmap_id >= 0)) { error_code |= - cvm::error("Error: mapName and mapID are mutually exclusive.\n"); + cvm::error("Error: mapName and mapID are mutually exclusive.\n", COLVARS_INPUT_ERROR); } // Parse optional group atoms = parse_group(conf, "atoms", true); - if (atoms != NULL) { + if (atoms) { // Using internal selection if (volmap_name.size()) { @@ -73,11 +57,11 @@ int colvar::map_total::init(std::string const &conf) if (volmap_id >= 0) { volmap_index = proxy->init_volmap_by_id(volmap_id); } - error_code |= volmap_index > 0 ? COLVARS_OK : COLVARS_INPUT_ERROR; + error_code |= (volmap_index >= 0) ? COLVARS_OK : COLVARS_INPUT_ERROR; } if (get_keyval(conf, "atomWeights", atom_weights, atom_weights)) { - if (atoms == NULL) { + if (!atoms) { error_code |= cvm::error("Error: weights can only be assigned when atoms " "are selected explicitly in Colvars.\n", COLVARS_INPUT_ERROR); @@ -132,11 +116,10 @@ void colvar::map_total::calc_gradients() void colvar::map_total::apply_force(colvarvalue const &force) { - colvarproxy *proxy = cvm::main()->proxy; if (atoms) { - if (!atoms->noforce) - atoms->apply_colvar_force(force.real_value); + cvc::apply_force(force); } else { + colvarproxy *proxy = cvm::main()->proxy; proxy->apply_volmap_force(volmap_index, force.real_value); } } diff --git a/src/colvardeps.h b/src/colvardeps.h index ac6a28851..1bd304b54 100644 --- a/src/colvardeps.h +++ b/src/colvardeps.h @@ -357,6 +357,8 @@ class colvardeps { f_cvc_lower_boundary, /// This CVC provides a default value for the colvar's upper boundary f_cvc_upper_boundary, + /// CVC accesses atom groups directly (as opposed to going throuh other objects) + f_cvc_explicit_atom_groups, /// CVC calculates atom gradients f_cvc_gradient, /// CVC calculates and stores explicit atom gradients on rank 0 diff --git a/src/colvarmodule_refs.h b/src/colvarmodule_refs.h index b3f584242..6fe6c65c7 100644 --- a/src/colvarmodule_refs.h +++ b/src/colvarmodule_refs.h @@ -441,14 +441,14 @@ feature_count_[std::string("orientationAngle colvar component (derived from orientation)")] = 0; feature_paper_map_[std::string("orientationAngle colvar component (derived from orientation)")] = "Fiorin2013"; - feature_count_[std::string("orientationProj colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("orientationProj colvar component (derived from orientation)")] = "Fiorin2013"; + feature_count_[std::string("orientationProj colvar component (derived from orientationAngle)")] = 0; + feature_paper_map_[std::string("orientationProj colvar component (derived from orientationAngle)")] = "Fiorin2013"; - feature_count_[std::string("spinAngle colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("spinAngle colvar component (derived from orientation)")] = "Fiorin2013"; + feature_count_[std::string("spinAngle colvar component (derived from tilt)")] = 0; + feature_paper_map_[std::string("spinAngle colvar component (derived from tilt)")] = "Fiorin2013"; - feature_count_[std::string("tilt colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("tilt colvar component (derived from orientation)")] = "Fiorin2013"; + feature_count_[std::string("tilt colvar component (derived from orientationProj)")] = 0; + feature_paper_map_[std::string("tilt colvar component (derived from orientationProj)")] = "Fiorin2013"; feature_count_[std::string("alpha colvar component")] = 0; feature_paper_map_[std::string("alpha colvar component")] = "Fiorin2013"; @@ -498,14 +498,14 @@ feature_count_[std::string("polarPhi colvar component")] = 0; feature_paper_map_[std::string("polarPhi colvar component")] = "Fu2017"; - feature_count_[std::string("eulerPhi colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("eulerPhi colvar component (derived from orientation)")] = "Fu2017"; + feature_count_[std::string("eulerPhi colvar component (derived from orientation_angle)")] = 0; + feature_paper_map_[std::string("eulerPhi colvar component (derived from orientation_angle)")] = "Fu2017"; - feature_count_[std::string("eulerTheta colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("eulerTheta colvar component (derived from orientation)")] = "Fu2017"; + feature_count_[std::string("eulerTheta colvar component (derived from orientation_angle)")] = 0; + feature_paper_map_[std::string("eulerTheta colvar component (derived from orientation_angle)")] = "Fu2017"; - feature_count_[std::string("eulerPsi colvar component (derived from orientation)")] = 0; - feature_paper_map_[std::string("eulerPsi colvar component (derived from orientation)")] = "Fu2017"; + feature_count_[std::string("eulerPsi colvar component (derived from orientation_angle)")] = 0; + feature_paper_map_[std::string("eulerPsi colvar component (derived from orientation_angle)")] = "Fu2017"; feature_count_[std::string("dipoleAngle colvar component")] = 0; feature_paper_map_[std::string("dipoleAngle colvar component")] = "Garate2019";