Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "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 : \verbatim
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 : \endverbatim
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 2527 : PLUMED_REGISTER_ACTION(MoleculeOrientation,"MOLECULES")
65 :
66 5 : void MoleculeOrientation::registerKeywords( Keywords& keys ) {
67 5 : VectorMultiColvar::registerKeywords( keys ); keys.use("VMEAN");
68 : 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 5 : "to be at the center of the vector connecting the first two atoms.");
72 5 : keys.reset_style("MOL","atoms");
73 5 : }
74 :
75 4 : MoleculeOrientation::MoleculeOrientation( const ActionOptions& ao ):
76 : Action(ao),
77 4 : VectorMultiColvar(ao)
78 : {
79 4 : int natoms=-1; std::vector<AtomNumber> all_atoms;
80 4 : readAtomsLikeKeyword("MOL",natoms,all_atoms);
81 4 : nvectors = std::floor( natoms / 2 );
82 4 : if( natoms%2!=0 && 2*nvectors+1!=natoms ) 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( natoms==2*nvectors+1 ) {
88 4 : std::vector<bool> catom_ind(natoms, false); catom_ind[natoms-1]=true;
89 4 : setAtomsForCentralAtom( catom_ind );
90 4 : }
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 20868 : void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
100 41733 : for(unsigned i=0; i<nvectors; ++i) {
101 20863 : Vector distance; distance=getSeparation( myatoms.getPosition(2*i+0), myatoms.getPosition(2*i+1) );
102 :
103 20865 : addAtomDerivatives( 2+3*i+0, 2*i+0, Vector(-1.0,0,0), myatoms );
104 20860 : addAtomDerivatives( 2+3*i+0, 2*i+1, Vector(+1.0,0,0), myatoms );
105 20862 : myatoms.addBoxDerivatives( 2+3*i+0, Tensor(distance,Vector(-1.0,0,0)) );
106 20863 : myatoms.addValue( 2+3*i+0, distance[0] );
107 :
108 20860 : addAtomDerivatives( 2+3*i+1, 2*i+0, Vector(0,-1.0,0), myatoms );
109 20870 : addAtomDerivatives( 2+3*i+1, 2*i+1, Vector(0,+1.0,0), myatoms );
110 20868 : myatoms.addBoxDerivatives( 2+3*i+1, Tensor(distance,Vector(0,-1.0,0)) );
111 20868 : myatoms.addValue( 2+3*i+1, distance[1] );
112 :
113 20865 : addAtomDerivatives( 2+3*i+2, 2*i+0, Vector(0,0,-1.0), myatoms );
114 20868 : addAtomDerivatives( 2+3*i+2, 2*i+1, Vector(0,0,+1.0), myatoms );
115 20869 : myatoms.addBoxDerivatives( 2+3*i+2, Tensor(distance,Vector(0,0,-1.0)) );
116 20866 : myatoms.addValue( 2+3*i+2, distance[2] );
117 : }
118 20870 : }
119 :
120 1041228 : void MoleculeOrientation::normalizeVector( std::vector<double>& vals ) const {
121 2082456 : for(unsigned i=0; i<nvectors; ++i) {
122 1041228 : double norm=0;
123 1041228 : 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 1041228 : 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 9984 : std::vector<double> weight( nvectors ), wdf( nvectors );
133 9984 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
134 4992 : double v=0; for(unsigned jcomp=0; jcomp<3; ++jcomp) v += myvals.get( 2+3*ivec+jcomp )*myvals.get( 2+3*ivec+jcomp );
135 4992 : v=sqrt(v); weight[ivec]=1.0; wdf[ivec]=1.0;
136 4992 : if( v>epsilon ) { weight[ivec] = 1.0 / v; wdf[ivec] = 1.0 / ( v*v*v ); }
137 : }
138 :
139 54672 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
140 49680 : unsigned jder=myvals.getActiveIndex(j);
141 99360 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
142 49680 : 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 198720 : for(unsigned jcomp=0; jcomp<3; ++jcomp) {
144 149040 : 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 4992 : }
148 4992 : }
149 :
150 : }
151 2523 : }
|
