Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "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 : 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 override;
59 : void calculateVector( multicolvar::AtomValuePack& myatoms ) const override;
60 : void normalizeVector( std::vector<double>& vals ) const override;
61 : void normalizeVectorDerivatives( MultiValue& myvals ) const override;
62 : };
63 :
64 13793 : PLUMED_REGISTER_ACTION(MoleculeOrientation,"MOLECULES")
65 :
66 8 : void MoleculeOrientation::registerKeywords( Keywords& keys ) {
67 8 : VectorMultiColvar::registerKeywords( keys );
68 8 : keys.use("MEAN");
69 8 : keys.use("VMEAN");
70 16 : keys.add("numbered","MOL","The numerical indices of the atoms in the molecule. The orientation of the molecule is equal to "
71 : "the vector connecting the first two atoms specified. If a third atom is specified its position "
72 : "is used to specify where the molecule is. If a third atom is not present the molecule is assumed "
73 : "to be at the center of the vector connecting the first two atoms.");
74 16 : keys.reset_style("MOL","atoms");
75 8 : }
76 :
77 4 : MoleculeOrientation::MoleculeOrientation( const ActionOptions& ao ):
78 : Action(ao),
79 4 : VectorMultiColvar(ao) {
80 : std::vector<AtomNumber> all_atoms;
81 8 : readAtomsLikeKeyword("MOL",-1,all_atoms);
82 4 : nvectors = std::floor( ablocks.size() / 2 );
83 4 : if( ablocks.size()%2!=0 && 2*nvectors+1!=ablocks.size() ) {
84 0 : error("number of atoms in molecule specification is wrong. Should be two or three.");
85 : }
86 :
87 4 : if( all_atoms.size()==0 ) {
88 0 : error("No atoms were specified");
89 : }
90 4 : setVectorDimensionality( 3*nvectors );
91 4 : setupMultiColvarBase( all_atoms );
92 :
93 4 : if( ablocks.size()==2*nvectors+1 ) {
94 4 : std::vector<bool> catom_ind(ablocks.size(), false);
95 4 : catom_ind[ablocks.size()-1]=true;
96 4 : setAtomsForCentralAtom( catom_ind );
97 : }
98 4 : }
99 :
100 0 : AtomNumber MoleculeOrientation::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
101 : plumed_dbg_assert( iatom<atom_lab.size() );
102 0 : plumed_assert( atom_lab[iatom].first==0 );
103 0 : return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
104 : }
105 :
106 16584 : void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
107 33168 : for(unsigned i=0; i<nvectors; ++i) {
108 16584 : Vector distance;
109 16584 : distance=getSeparation( myatoms.getPosition(2*i+0), myatoms.getPosition(2*i+1) );
110 :
111 16584 : addAtomDerivatives( 2+3*i+0, 2*i+0, Vector(-1.0,0,0), myatoms );
112 16584 : addAtomDerivatives( 2+3*i+0, 2*i+1, Vector(+1.0,0,0), myatoms );
113 16584 : myatoms.addBoxDerivatives( 2+3*i+0, Tensor(distance,Vector(-1.0,0,0)) );
114 16584 : myatoms.addValue( 2+3*i+0, distance[0] );
115 :
116 16584 : addAtomDerivatives( 2+3*i+1, 2*i+0, Vector(0,-1.0,0), myatoms );
117 16584 : addAtomDerivatives( 2+3*i+1, 2*i+1, Vector(0,+1.0,0), myatoms );
118 16584 : myatoms.addBoxDerivatives( 2+3*i+1, Tensor(distance,Vector(0,-1.0,0)) );
119 16584 : myatoms.addValue( 2+3*i+1, distance[1] );
120 :
121 16584 : addAtomDerivatives( 2+3*i+2, 2*i+0, Vector(0,0,-1.0), myatoms );
122 16584 : addAtomDerivatives( 2+3*i+2, 2*i+1, Vector(0,0,+1.0), myatoms );
123 16584 : myatoms.addBoxDerivatives( 2+3*i+2, Tensor(distance,Vector(0,0,-1.0)) );
124 16584 : myatoms.addValue( 2+3*i+2, distance[2] );
125 : }
126 16584 : }
127 :
128 1041228 : void MoleculeOrientation::normalizeVector( std::vector<double>& vals ) const {
129 2082456 : for(unsigned i=0; i<nvectors; ++i) {
130 : double norm=0;
131 4164912 : for(unsigned j=0; j<3; ++j) {
132 3123684 : norm += vals[2+3*i+j]*vals[2+3*i+j];
133 : }
134 1041228 : norm = sqrt(norm);
135 :
136 : double inorm = 1.0;
137 1041228 : if( norm>epsilon ) {
138 1041228 : inorm = 1.0 / norm;
139 : }
140 4164912 : for(unsigned j=0; j<3; ++j) {
141 3123684 : vals[2+3*i+j] = inorm*vals[2+3*i+j];
142 : }
143 : }
144 1041228 : }
145 :
146 4992 : void MoleculeOrientation::normalizeVectorDerivatives( MultiValue& myvals ) const {
147 4992 : std::vector<double> weight( nvectors ), wdf( nvectors );
148 9984 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
149 : double v=0;
150 19968 : for(unsigned jcomp=0; jcomp<3; ++jcomp) {
151 14976 : v += myvals.get( 2+3*ivec+jcomp )*myvals.get( 2+3*ivec+jcomp );
152 : }
153 4992 : v=sqrt(v);
154 4992 : weight[ivec]=1.0;
155 4992 : wdf[ivec]=1.0;
156 4992 : if( v>epsilon ) {
157 4992 : weight[ivec] = 1.0 / v;
158 4992 : wdf[ivec] = 1.0 / ( v*v*v );
159 : }
160 : }
161 :
162 54672 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
163 49680 : unsigned jder=myvals.getActiveIndex(j);
164 99360 : for(unsigned ivec=0; ivec<nvectors; ++ivec) {
165 : double comp2=0.0;
166 198720 : for(unsigned jcomp=0; jcomp<3; ++jcomp) {
167 149040 : comp2 += myvals.get(2+3*ivec+jcomp)*myvals.getDerivative( 2+3*ivec+jcomp, jder );
168 : }
169 198720 : for(unsigned jcomp=0; jcomp<3; ++jcomp) {
170 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) );
171 : }
172 : }
173 : }
174 4992 : }
175 :
176 : }
177 : }
|
