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() { return true; }
93 0 : void retrieveDomain( std::string& min, std::string& max ) { min="-pi"; max="+pi"; }
94 : };
95 :
96 10419 : PLUMED_REGISTER_ACTION(InterMolecularTorsions,"INTERMOLECULARTORSIONS")
97 :
98 1 : void InterMolecularTorsions::registerKeywords( Keywords& keys ) {
99 1 : MultiColvarBase::registerKeywords( keys );
100 2 : 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");
101 2 : 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. "
102 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
103 2 : 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. "
104 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
105 2 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
106 2 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
107 2 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
108 2 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
109 2 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
110 : "The following provides information on the \\ref switchingfunction that are available. "
111 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
112 : // Use actionWithDistributionKeywords
113 1 : keys.remove("LOWMEM");
114 2 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
115 1 : }
116 :
117 0 : InterMolecularTorsions::InterMolecularTorsions(const ActionOptions& ao):
118 : Action(ao),
119 0 : MultiColvarBase(ao)
120 : {
121 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
122 0 : if( getBaseMultiColvar(i)->getNumberOfQuantities()!=5 ) error("input multicolvar does not calculate molecular orientations");
123 : }
124 : // The weight of this does have derivatives
125 0 : weightHasDerivatives=true;
126 :
127 : // Read in the switching function
128 0 : std::string sw, errors; parse("SWITCH",sw);
129 0 : if(sw.length()>0) {
130 0 : switchingFunction.set(sw,errors);
131 : } else {
132 0 : double r_0=-1.0, d_0; int nn, mm;
133 0 : parse("NN",nn); parse("MM",mm);
134 0 : parse("R_0",r_0); parse("D_0",d_0);
135 0 : if( r_0<0.0 ) error("you must set a value for R_0");
136 0 : switchingFunction.set(nn,mm,r_0,d_0);
137 : }
138 0 : log.printf(" calculating number of links with atoms separation of %s\n",( switchingFunction.description() ).c_str() );
139 0 : std::vector<AtomNumber> all_atoms; readTwoGroups( "MOLS", "MOLSA", "MOLSB", all_atoms );
140 0 : setupMultiColvarBase( all_atoms ); setLinkCellCutoff( switchingFunction.get_dmax() );
141 :
142 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
143 0 : if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
144 : }
145 :
146 : // Create holders for the collective variable
147 0 : readVesselKeywords();
148 0 : plumed_assert( getNumberOfVessels()==0 );
149 0 : std::string input; addVessel( "SUM", input, -1 );
150 0 : readVesselKeywords();
151 0 : }
152 :
153 0 : double InterMolecularTorsions::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
154 0 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
155 0 : double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
156 :
157 0 : if( !doNotCalculateDerivatives() ) {
158 0 : addAtomDerivatives( 0, 0, (-dfunc)*weight*distance, myatoms );
159 0 : addAtomDerivatives( 0, 1, (dfunc)*weight*distance, myatoms );
160 0 : myatoms.addBoxDerivatives( 0, (-dfunc)*weight*Tensor(distance,distance) );
161 : }
162 0 : return sw;
163 : }
164 :
165 0 : double InterMolecularTorsions::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
166 0 : Vector v1, v2, dv1, dv2, dconn, conn = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
167 :
168 : // Retrieve vectors
169 0 : std::vector<double> orient0( 5 ), orient1( 5 );
170 0 : getInputData( 0, true, myatoms, orient0 );
171 0 : getInputData( 1, true, myatoms, orient1 );
172 0 : for(unsigned i=0; i<3; ++i) { v1[i]=orient0[2+i]; v2[i]=orient1[2+i]; }
173 0 : if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) return 1.0;
174 :
175 : // Evaluate angle
176 0 : Torsion t; double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
177 0 : for(unsigned i=0; i<3; ++i) { orient0[i+2]=dv1[i]; orient1[i+2]=dv2[i]; }
178 :
179 : // And accumulate derivatives
180 0 : if( !doNotCalculateDerivatives() ) {
181 0 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
182 0 : mergeInputDerivatives( 1, 2, orient1.size(), 0, orient0, myder0, myatoms );
183 0 : MultiValue& myder1=getInputDerivatives( 1, true, myatoms );
184 0 : mergeInputDerivatives( 1, 2, orient0.size(), 1, orient1, myder1, myatoms );
185 0 : addAtomDerivatives( 1, 0, -dconn, myatoms ); addAtomDerivatives( 1, 1, dconn, myatoms );
186 0 : myatoms.addBoxDerivatives( 1, -extProduct( conn, dconn ) );
187 : }
188 :
189 : return angle;
190 : }
191 :
192 : }
193 : }
|
