Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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/MultiColvarFunction.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/SwitchingFunction.h"
25 : #include "tools/Torsion.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MCOLVARF INTERMOLECULARTORSIONS
31 : /*
32 : Calculate torsions between axis of adjacent molecules
33 :
34 : This command can be used to calculate the intermolecular torsional angles between the orientations of nearby molecules. The orientation of a
35 : molecule can be calculated by using either the \ref MOLECULES or the \ref PLANES commands. These two commands calculate the orientation of a
36 : 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
37 : molecules as objects that lie at a point in space and that have an orientation. This command calculates the torsional angles between the orientations
38 : 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
39 : 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
40 : 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
41 : that are close together.
42 :
43 : \par Examples
44 :
45 : The example input below is necessarily but gives you an idea of what can be achieved using this action.
46 : The orientations and positions of four molecules are defined using the \ref MOLECULES action as the position of the
47 : centeres of mass of the two atoms specified and the direction of the vector connecting the two atoms that were specified.
48 : The torsional angles between the molecules are then calculated by the \ref INTERMOLECULARTORSIONS command labelled tt_p.
49 : We then compute a \ref HISTOGRAM that shows the distribution that these torsional angles take in the structure. The weight
50 : a given torsional angle contributes to this \ref HISTOGRAM is determined using a \ref switchingfunction that acts on the distance
51 : between the two molecules. As such the torsional angles between molecules that are close together contribute a high weight to the
52 : histogram while the torsional angles between molecules that are far apart does not contribute to the histogram. The histogram is
53 : averaged over the whole trajectory and output once all the trajectory frames have been read.
54 :
55 : \verbatim
56 : m1: MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8
57 : tt_p: INTERMOLECULARTORSIONS MOLS=m1 SWITCH={RATIONAL R_0=0.25 D_0=2.0 D_MAX=3.0}
58 : htt_p: HISTOGRAM DATA=tt_p GRID_MIN=-pi GRID_MAX=pi BANDWIDTH=0.1 GRID_BIN=200 STRIDE=1
59 : DUMPGRID GRID=htt_p FILE=myhist.out
60 : \endverbatim
61 :
62 : */
63 : //+ENDPLUMEDOC
64 :
65 : namespace PLMD {
66 : namespace crystallization {
67 :
68 0 : class InterMolecularTorsions : public multicolvar::MultiColvarFunction {
69 : private:
70 : /// The switching function that tells us if atoms are close enough together
71 : SwitchingFunction switchingFunction;
72 : public:
73 : static void registerKeywords( Keywords& keys );
74 : explicit InterMolecularTorsions(const ActionOptions&);
75 : /// Do the stuff with the switching functions
76 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const ;
77 : /// Actually do the calculation
78 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
79 : /// Is the variable periodic
80 0 : bool isPeriodic() { return true; }
81 0 : void retrieveDomain( std::string& min, std::string& max ) { min="-pi"; max="+pi"; }
82 : };
83 :
84 2523 : PLUMED_REGISTER_ACTION(InterMolecularTorsions,"INTERMOLECULARTORSIONS")
85 :
86 1 : void InterMolecularTorsions::registerKeywords( Keywords& keys ) {
87 1 : MultiColvarFunction::registerKeywords( keys );
88 1 : 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");
89 : keys.add("atoms-1","MOLSA","In this version of the input the torsional angles between all pairs of atoms including one atom from MOLA one atom from MOLB will be computed. "
90 1 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
91 : keys.add("atoms-1","MOLSB","In this version of the input the torsional angles between all pairs of atoms including one atom from MOLA one atom from MOLB will be computed. "
92 1 : "This should be the label/s of \\ref MOLECULES or \\ref PLANES actions");
93 1 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
94 1 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
95 1 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
96 1 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
97 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
98 : "The following provides information on the \\ref switchingfunction that are available. "
99 1 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
100 : // Use actionWithDistributionKeywords
101 1 : keys.remove("LOWMEM"); keys.remove("DATA");
102 1 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
103 1 : }
104 :
105 0 : InterMolecularTorsions::InterMolecularTorsions(const ActionOptions& ao):
106 : Action(ao),
107 0 : MultiColvarFunction(ao)
108 : {
109 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
110 0 : if( getBaseMultiColvar(i)->getNumberOfQuantities()!=5 ) error("input multicolvar does not calculate molecular orientations");
111 : }
112 : // The weight of this does have derivatives
113 0 : weightHasDerivatives=true;
114 :
115 : // Read in the switching function
116 0 : std::string sw, errors; parse("SWITCH",sw);
117 0 : if(sw.length()>0) {
118 0 : switchingFunction.set(sw,errors);
119 : } else {
120 0 : double r_0=-1.0, d_0; int nn, mm;
121 0 : parse("NN",nn); parse("MM",mm);
122 0 : parse("R_0",r_0); parse("D_0",d_0);
123 0 : if( r_0<0.0 ) error("you must set a value for R_0");
124 0 : switchingFunction.set(nn,mm,r_0,d_0);
125 : }
126 0 : log.printf(" calculating number of links with atoms separation of %s\n",( switchingFunction.description() ).c_str() );
127 0 : std::vector<AtomNumber> all_atoms; readTwoGroups( "MOLS", "MOLSA", "MOLSB", all_atoms );
128 0 : setupMultiColvarBase( all_atoms ); setLinkCellCutoff( switchingFunction.get_dmax() );
129 :
130 0 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
131 0 : if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
132 : }
133 :
134 : // Create holders for the collective variable
135 0 : readVesselKeywords();
136 0 : plumed_assert( getNumberOfVessels()==0 );
137 0 : std::string input; addVessel( "SUM", input, -1 );
138 0 : readVesselKeywords();
139 0 : }
140 :
141 0 : double InterMolecularTorsions::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
142 0 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
143 0 : double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
144 :
145 0 : if( !doNotCalculateDerivatives() ) {
146 0 : addAtomDerivatives( 0, 0, (-dfunc)*weight*distance, myatoms );
147 0 : addAtomDerivatives( 0, 1, (dfunc)*weight*distance, myatoms );
148 0 : myatoms.addBoxDerivatives( 0, (-dfunc)*weight*Tensor(distance,distance) );
149 : }
150 0 : return sw;
151 : }
152 :
153 0 : double InterMolecularTorsions::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
154 0 : Vector v1, v2, dv1, dv2, dconn, conn = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
155 :
156 : // Retrieve vectors
157 0 : std::vector<double> orient0( 5 ), orient1( 5 );
158 0 : getInputData( 0, true, myatoms, orient0 );
159 0 : getInputData( 1, true, myatoms, orient1 );
160 0 : for(unsigned i=0; i<3; ++i) { v1[i]=orient0[2+i]; v2[i]=orient1[2+i]; }
161 0 : if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) return 1.0;
162 :
163 : // Evaluate angle
164 0 : Torsion t; double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 );
165 0 : for(unsigned i=0; i<3; ++i) { orient0[i+2]=dv1[i]; orient1[i+2]=dv2[i]; }
166 :
167 : // And accumulate derivatives
168 0 : if( !doNotCalculateDerivatives() ) {
169 0 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
170 0 : mergeInputDerivatives( 1, 2, orient1.size(), 0, orient0, myder0, myatoms );
171 0 : MultiValue& myder1=getInputDerivatives( 1, true, myatoms );
172 0 : mergeInputDerivatives( 1, 2, orient0.size(), 1, orient1, myder1, myatoms );
173 0 : addAtomDerivatives( 1, 0, -dconn, myatoms ); addAtomDerivatives( 1, 1, dconn, myatoms );
174 0 : myatoms.addBoxDerivatives( 1, -extProduct( conn, dconn ) );
175 : }
176 :
177 0 : return angle;
178 : }
179 :
180 : }
181 2523 : }
|
