Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 "core/ActionRegister.h"
23 : #include "VectorMultiColvar.h"
24 :
25 : namespace PLMD {
26 : namespace crystallization {
27 :
28 : //+PLUMEDOC MCOLVAR MOLECULES
29 : /*
30 : Calculate the vectors connecting a pair of atoms in order to represent the orientation of a molecule.
31 :
32 : At its simplest this command can be used to calculate the average length of an internal vector in a
33 : collection of different molecules. When used in conjunction with MutiColvarFunctions in can be used
34 : to do a variety of more complex tasks.
35 :
36 : \par Examples
37 :
38 : The following input tells plumed to calculate the distances between two of the atoms in a molecule.
39 : This is done for the same set of atoms four different molecules and the average separation is then
40 : calculated.
41 :
42 : \plumedfile
43 : MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8 MEAN LABEL=mm
44 : PRINT ARG=mm.mean FILE=colvar
45 : \endplumedfile
46 :
47 :
48 : */
49 : //+ENDPLUMEDOC
50 :
51 :
52 8 : class MoleculeOrientation : public VectorMultiColvar {
53 : private:
54 : unsigned nvectors;
55 : public:
56 : static void registerKeywords( Keywords& keys );
57 : explicit MoleculeOrientation( const ActionOptions& ao );
58 : AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const ;
59 : void calculateVector( multicolvar::AtomValuePack& myatoms ) const;
60 : void normalizeVector( std::vector<double>& vals ) const ;
61 : void normalizeVectorDerivatives( MultiValue& myvals ) const ;
62 : };
63 :
64 6456 : PLUMED_REGISTER_ACTION(MoleculeOrientation,"MOLECULES")
65 :
66 5 : void MoleculeOrientation::registerKeywords( Keywords& keys ) {
67 10 : VectorMultiColvar::registerKeywords( keys ); keys.use("VMEAN");
68 20 : keys.add("numbered","MOL","The numerical indices of the atoms in the molecule. The orientation of the molecule is equal to "
69 : "the vector connecting the first two atoms specified. If a third atom is specified its position "
70 : "is used to specify where the molecule is. If a third atom is not present the molecule is assumed "
71 : "to be at the center of the vector connecting the first two atoms.");
72 15 : keys.reset_style("MOL","atoms");
73 5 : }
74 :
75 4 : MoleculeOrientation::MoleculeOrientation( const ActionOptions& ao ):
76 : Action(ao),
77 4 : VectorMultiColvar(ao)
78 : {
79 : std::vector<AtomNumber> all_atoms;
80 8 : readAtomsLikeKeyword("MOL",-1,all_atoms);
81 8 : nvectors = std::floor( ablocks.size() / 2 );
82 4 : if( ablocks.size()%2!=0 && 2*nvectors+1!=ablocks.size() ) error("number of atoms in molecule specification is wrong. Should be two or three.");
83 :
84 4 : if( all_atoms.size()==0 ) error("No atoms were specified");
85 4 : setVectorDimensionality( 3*nvectors ); setupMultiColvarBase( all_atoms );
86 :
87 4 : if( ablocks.size()==2*nvectors+1 ) {
88 8 : std::vector<bool> catom_ind(ablocks.size(), false); catom_ind[ablocks.size()-1]=true;
89 4 : setAtomsForCentralAtom( catom_ind );
90 : }
91 4 : }
92 :
93 0 : AtomNumber MoleculeOrientation::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
94 : plumed_dbg_assert( iatom<atom_lab.size() );
95 0 : plumed_assert( atom_lab[iatom].first==0 );
96 0 : return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
97 : }
98 :
99 16584 : void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
100 49752 : for(unsigned i=0; i<nvectors; ++i) {
101 33168 : Vector distance; distance=getSeparation( myatoms.getPosition(2*i+0), myatoms.getPosition(2*i+1) );
102 :
103 16584 : addAtomDerivatives( 2+3*i+0, 2*i+0, Vector(-1.0,0,0), myatoms );
104 16584 : addAtomDerivatives( 2+3*i+0, 2*i+1, Vector(+1.0,0,0), myatoms );
105 16584 : myatoms.addBoxDerivatives( 2+3*i+0, Tensor(distance,Vector(-1.0,0,0)) );
106 16584 : myatoms.addValue( 2+3*i+0, distance[0] );
107 :
108 16584 : addAtomDerivatives( 2+3*i+1, 2*i+0, Vector(0,-1.0,0), myatoms );
109 16584 : addAtomDerivatives( 2+3*i+1, 2*i+1, Vector(0,+1.0,0), myatoms );
110 16584 : myatoms.addBoxDerivatives( 2+3*i+1, Tensor(distance,Vector(0,-1.0,0)) );
111 16584 : myatoms.addValue( 2+3*i+1, distance[1] );
112 :
113 16584 : addAtomDerivatives( 2+3*i+2, 2*i+0, Vector(0,0,-1.0), myatoms );
114 16584 : addAtomDerivatives( 2+3*i+2, 2*i+1, Vector(0,0,+1.0), myatoms );
115 16584 : myatoms.addBoxDerivatives( 2+3*i+2, Tensor(distance,Vector(0,0,-1.0)) );
116 16584 : myatoms.addValue( 2+3*i+2, distance[2] );
117 : }
118 16584 : }
119 :
120 1041228 : void MoleculeOrientation::normalizeVector( std::vector<double>& vals ) const {
121 3123684 : for(unsigned i=0; i<nvectors; ++i) {
122 : double norm=0;
123 7288596 : for(unsigned j=0; j<3; ++j) norm += vals[2+3*i+j]*vals[2+3*i+j];
124 1041228 : norm = sqrt(norm);
125 :
126 1041228 : double inorm = 1.0; if( norm>epsilon ) inorm = 1.0 / norm;
127 7288596 : for(unsigned j=0; j<3; ++j) vals[2+3*i+j] = inorm*vals[2+3*i+j];
128 : }
129 1041228 : }
130 :
131 4992 : void MoleculeOrientation::normalizeVectorDerivatives( MultiValue& myvals ) const {
132 4992 : std::vector<double> weight( nvectors ), wdf( nvectors );
133 14976 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
134 34944 : double v=0; for(unsigned jcomp=0; jcomp<3; ++jcomp) v += myvals.get( 2+3*ivec+jcomp )*myvals.get( 2+3*ivec+jcomp );
135 14976 : v=sqrt(v); weight[ivec]=1.0; wdf[ivec]=1.0;
136 14976 : if( v>epsilon ) { weight[ivec] = 1.0 / v; wdf[ivec] = 1.0 / ( v*v*v ); }
137 : }
138 :
139 104352 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
140 : unsigned jder=myvals.getActiveIndex(j);
141 149040 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
142 347760 : double comp2=0.0; for(unsigned jcomp=0; jcomp<3; ++jcomp) comp2 += myvals.get(2+3*ivec+jcomp)*myvals.getDerivative( 2+3*ivec+jcomp, jder );
143 347760 : for(unsigned jcomp=0; jcomp<3; ++jcomp) {
144 745200 : myvals.setDerivative( 2+3*ivec+jcomp, jder, weight[ivec]*myvals.getDerivative( 2+3*ivec+jcomp, jder ) - wdf[ivec]*comp2*myvals.get(2+3*ivec+jcomp) );
145 : }
146 : }
147 : }
148 4992 : }
149 :
150 : }
151 4839 : }
|
