Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "multicolvar/MultiColvarBase.h"
23 : #include "multicolvar/AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "tools/Torsion.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : //+PLUMEDOC MCOLVARF INTERMOLECULARTORSIONS
32 : /*
33 : Calculate torsion angles between vectors on adjacent molecules
34 :
35 : This variable can be used to calculate the average torsional angles between vectors. In other words,
36 : it can be used to compute quantities like this:
37 :
38 : \f[
39 : s = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \theta_{ij} }{ \sum_{i \ne j} \sigma(r_{ij}) }
40 : \f]
41 :
42 : Here the sums run over all pairs of molecules. \f$\sigma(r_{ij})\f$ is a \ref switchingfunction that
43 : action on the distance between the centers of molecules \f$i\f$ and \f$j\f$. \f$\theta_{ij}\f$ is then
44 : the torsional angle between an orientation vector for molecule \f$i\f$ and molecule \f$j\f$.
45 :
46 : This command can be used to calculate the intermolecular torsional angles between the orientations of nearby molecules. The orientation of a
47 : molecule can be calculated by using either the \ref MOLECULES or the \ref PLANES commands. These two commands calculate the orientation of a
48 : bond in the molecule or the orientation of a plane containing three of the molecule's atoms. Furthermore, when we use these commands we think of
49 : molecules as objects that lie at a point in space and that have an orientation. This command calculates the torsional angles between the orientations
50 : of these objects. We can then calculates functions of a large number of these torsional angles that measures things such as the number of torsional
51 : angles that are within a particular range. Because it is often useful to only consider the torsional angles between objects that are within a certain
52 : distance of each other we can, when calculating these sums, perform a weighted sum and use a \ref switchingfunction to ensure that we focus on molecules
53 : that are close together.
54 :
55 : \par Examples
56 :
57 : The example input below is necessarily but gives you an idea of what can be achieved using this action.
58 : The orientations and positions of four molecules are defined using the \ref MOLECULES action as the position of the
59 : centers of mass of the two atoms specified and the direction of the vector connecting the two atoms that were specified.
60 : The torsional angles between the molecules are then calculated by the \ref INTERMOLECULARTORSIONS command labelled torsion_p.
61 : We then compute a \ref HISTOGRAM that shows the distribution that these torsional angles take in the structure. The weight
62 : a given torsional angle contributes to this \ref HISTOGRAM is determined using a \ref switchingfunction that acts on the distance
63 : between the two molecules. As such the torsional angles between molecules that are close together contribute a high weight to the
64 : histogram while the torsional angles between molecules that are far apart does not contribute to the histogram. The histogram is
65 : averaged over the whole trajectory and output once all the trajectory frames have been read.
66 :
67 : \plumedfile
68 : m1: MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8
69 : torsion_p: INTERMOLECULARTORSIONS MOLS=m1 SWITCH={RATIONAL R_0=0.25 D_0=2.0 D_MAX=3.0}
70 : htt_p: HISTOGRAM DATA=torsion_p GRID_MIN=-pi GRID_MAX=pi BANDWIDTH=0.1 GRID_BIN=200 STRIDE=1
71 : DUMPGRID GRID=htt_p FILE=myhist.out
72 : \endplumedfile
73 :
74 : */
75 : //+ENDPLUMEDOC
76 :
77 : namespace PLMD {
78 : namespace crystallization {
79 :
80 : class InterMolecularTorsions : public multicolvar::MultiColvarBase {
81 : private:
82 : /// The switching function that tells us if atoms are close enough together
83 : SwitchingFunction switchingFunction;
84 : public:
85 : static void registerKeywords( Keywords& keys );
86 : explicit InterMolecularTorsions(const ActionOptions&);
87 : /// Do the stuff with the switching functions
88 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const ;
89 : /// Actually do the calculation
90 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
91 : /// Is the variable periodic
92 0 : bool isPeriodic() {
93 0 : return true;
94 : }
95 0 : void retrieveDomain( std::string& min, std::string& max ) {
96 : min="-pi";
97 : max="+pi";
98 0 : }
99 : };
100 :
101 13785 : PLUMED_REGISTER_ACTION(InterMolecularTorsions,"INTERMOLECULARTORSIONS")
102 :
103 4 : void InterMolecularTorsions::registerKeywords( Keywords& keys ) {
104 4 : MultiColvarBase::registerKeywords( keys );
105 8 : keys.add("atoms","MOLS","The molecules you would like to calculate the torsional angles between. This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
106 8 : keys.add("atoms-1","MOLSA","In this version of the input the torsional angles between all pairs of atoms including one atom from MOLSA one atom from MOLSB will be computed. "
107 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
108 8 : keys.add("atoms-1","MOLSB","In this version of the input the torsional angles between all pairs of atoms including one atom from MOLSA one atom from MOLSB will be computed. "
109 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
110 8 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
111 8 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
112 8 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
113 8 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
114 8 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
115 : "The following provides information on the \\ref switchingfunction that are available. "
116 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
117 : // Use actionWithDistributionKeywords
118 4 : keys.remove("LOWMEM");
119 8 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
120 4 : }
121 :
122 0 : InterMolecularTorsions::InterMolecularTorsions(const ActionOptions& ao):
123 : Action(ao),
124 0 : MultiColvarBase(ao) {
125 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
126 0 : if( getBaseMultiColvar(i)->getNumberOfQuantities()!=5 ) {
127 0 : error("input multicolvar does not calculate molecular orientations");
128 : }
129 : }
130 : // The weight of this does have derivatives
131 0 : weightHasDerivatives=true;
132 :
133 : // Read in the switching function
134 : std::string sw, errors;
135 0 : parse("SWITCH",sw);
136 0 : if(sw.length()>0) {
137 0 : switchingFunction.set(sw,errors);
138 : } else {
139 0 : double r_0=-1.0, d_0;
140 : int nn, mm;
141 0 : parse("NN",nn);
142 0 : parse("MM",mm);
143 0 : parse("R_0",r_0);
144 0 : parse("D_0",d_0);
145 0 : if( r_0<0.0 ) {
146 0 : error("you must set a value for R_0");
147 : }
148 0 : switchingFunction.set(nn,mm,r_0,d_0);
149 : }
150 0 : log.printf(" calculating number of links with atoms separation of %s\n",( switchingFunction.description() ).c_str() );
151 : std::vector<AtomNumber> all_atoms;
152 0 : readTwoGroups( "MOLS", "MOLSA", "MOLSB", all_atoms );
153 0 : setupMultiColvarBase( all_atoms );
154 0 : setLinkCellCutoff( switchingFunction.get_dmax() );
155 :
156 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
157 0 : if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) {
158 0 : error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
159 : }
160 : }
161 :
162 : // Create holders for the collective variable
163 0 : readVesselKeywords();
164 0 : plumed_assert( getNumberOfVessels()==0 );
165 : std::string input;
166 0 : addVessel( "SUM", input, -1 );
167 0 : readVesselKeywords();
168 0 : }
169 :
170 0 : double InterMolecularTorsions::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
171 0 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
172 0 : double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
173 :
174 0 : if( !doNotCalculateDerivatives() ) {
175 0 : addAtomDerivatives( 0, 0, (-dfunc)*weight*distance, myatoms );
176 0 : addAtomDerivatives( 0, 1, (dfunc)*weight*distance, myatoms );
177 0 : myatoms.addBoxDerivatives( 0, (-dfunc)*weight*Tensor(distance,distance) );
178 : }
179 0 : return sw;
180 : }
181 :
182 0 : double InterMolecularTorsions::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
183 0 : Vector v1, v2, dv1, dv2, dconn, conn = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
184 :
185 : // Retrieve vectors
186 0 : std::vector<double> orient0( 5 ), orient1( 5 );
187 0 : getInputData( 0, true, myatoms, orient0 );
188 0 : getInputData( 1, true, myatoms, orient1 );
189 0 : for(unsigned i=0; i<3; ++i) {
190 0 : v1[i]=orient0[2+i];
191 0 : v2[i]=orient1[2+i];
192 : }
193 0 : if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) {
194 : return 1.0;
195 : }
196 :
197 : // Evaluate angle
198 : Torsion t;
199 0 : double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
200 0 : for(unsigned i=0; i<3; ++i) {
201 0 : orient0[i+2]=dv1[i];
202 0 : orient1[i+2]=dv2[i];
203 : }
204 :
205 : // And accumulate derivatives
206 0 : if( !doNotCalculateDerivatives() ) {
207 0 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
208 0 : mergeInputDerivatives( 1, 2, orient1.size(), 0, orient0, myder0, myatoms );
209 0 : MultiValue& myder1=getInputDerivatives( 1, true, myatoms );
210 0 : mergeInputDerivatives( 1, 2, orient0.size(), 1, orient1, myder1, myatoms );
211 0 : addAtomDerivatives( 1, 0, -dconn, myatoms );
212 0 : addAtomDerivatives( 1, 1, dconn, myatoms );
213 0 : myatoms.addBoxDerivatives( 1, -extProduct( conn, dconn ) );
214 : }
215 :
216 : return angle;
217 : }
218 :
219 : }
220 : }
|
