Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 PLANES
29 : /*
30 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule.
31 :
32 : \par Examples
33 :
34 :
35 : */
36 : //+ENDPLUMEDOC
37 :
38 :
39 : class MoleculePlane : public VectorMultiColvar {
40 : private:
41 : public:
42 : static void registerKeywords( Keywords& keys );
43 : explicit MoleculePlane( const ActionOptions& ao );
44 : AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const override;
45 : void calculateVector( multicolvar::AtomValuePack& myatoms ) const override;
46 : };
47 :
48 13789 : PLUMED_REGISTER_ACTION(MoleculePlane,"PLANES")
49 :
50 6 : void MoleculePlane::registerKeywords( Keywords& keys ) {
51 6 : VectorMultiColvar::registerKeywords( keys );
52 6 : keys.use("VMEAN");
53 12 : keys.add("numbered","MOL","The numerical indices of the atoms in the molecule. If three atoms are specified the orientation "
54 : "of the molecule is taken as the normal to the plane containing the vector connecting the first and "
55 : "second atoms and the vector connecting the second and third atoms. If four atoms are specified the "
56 : "orientation of the molecule is taken as the normal to the plane containing the vector connecting the "
57 : "first and second atoms and the vector connecting the third and fourth atoms. The molecule is always "
58 : "assumed to lie at the geometric center for the three/four atoms.");
59 12 : keys.reset_style("MOL","atoms");
60 6 : }
61 :
62 2 : MoleculePlane::MoleculePlane( const ActionOptions& ao ):
63 : Action(ao),
64 2 : VectorMultiColvar(ao) {
65 : std::vector<AtomNumber> all_atoms;
66 4 : readAtomsLikeKeyword("MOL",-1,all_atoms);
67 2 : if( ablocks.size()!=3 && ablocks.size()!=4 ) {
68 0 : error("number of atoms in molecule specification is wrong. Should be three or four.");
69 : }
70 :
71 2 : if( all_atoms.size()==0 ) {
72 0 : error("No atoms were specified");
73 : }
74 2 : setVectorDimensionality( 3 );
75 2 : setupMultiColvarBase( all_atoms );
76 2 : }
77 :
78 0 : AtomNumber MoleculePlane::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
79 : plumed_dbg_assert( iatom<atom_lab.size() );
80 0 : plumed_assert( atom_lab[iatom].first==0 );
81 0 : return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
82 : }
83 :
84 8888 : void MoleculePlane::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
85 8888 : Vector d1, d2, cp;
86 8888 : if( myatoms.getNumberOfAtoms()==3 ) {
87 8888 : d1=getSeparation( myatoms.getPosition(1), myatoms.getPosition(0) );
88 8888 : d2=getSeparation( myatoms.getPosition(1), myatoms.getPosition(2) );
89 : } else {
90 0 : d1=getSeparation( myatoms.getPosition(1), myatoms.getPosition(0) );
91 0 : d2=getSeparation( myatoms.getPosition(2), myatoms.getPosition(3) );
92 : }
93 8888 : cp = crossProduct( d1, d2 );
94 :
95 8888 : addAtomDerivatives( 2, 0, crossProduct( Vector(-1.0,0,0), d2 ), myatoms );
96 8888 : if( myatoms.getNumberOfAtoms()==3 ) {
97 8888 : addAtomDerivatives( 2, 1, crossProduct( Vector(+1.0,0,0), d2 ) + crossProduct( Vector(-1.0,0,0), d1 ), myatoms );
98 8888 : addAtomDerivatives( 2, 2, crossProduct( Vector(+1.0,0,0), d1 ), myatoms );
99 : } else {
100 0 : addAtomDerivatives( 2, 1, crossProduct( Vector(+1.0,0,0), d2 ), myatoms );
101 0 : addAtomDerivatives( 2, 2, crossProduct( Vector(-1.0,0,0), d1 ), myatoms );
102 0 : addAtomDerivatives( 2, 3, crossProduct( Vector(+1.0,0,0), d1 ), myatoms );
103 : }
104 8888 : myatoms.addBoxDerivatives( 2, Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1)) );
105 8888 : myatoms.addValue( 2, cp[0] );
106 :
107 8888 : addAtomDerivatives( 3, 0, crossProduct( Vector(0,-1.0,0), d2 ), myatoms );
108 8888 : if( myatoms.getNumberOfAtoms()==3 ) {
109 8888 : addAtomDerivatives( 3, 1, crossProduct( Vector(0,+1.0,0), d2 ) + crossProduct( Vector(0,-1.0,0), d1 ), myatoms );
110 8888 : addAtomDerivatives( 3, 2, crossProduct( Vector(0,+1.0,0), d1 ), myatoms );
111 : } else {
112 0 : addAtomDerivatives( 3, 1, crossProduct( Vector(0,+1.0,0), d2 ), myatoms );
113 0 : addAtomDerivatives( 3, 2, crossProduct( Vector(0,-1.0,0), d1 ), myatoms );
114 0 : addAtomDerivatives( 3, 3, crossProduct( Vector(0,+1.0,0), d1 ), myatoms );
115 : }
116 8888 : myatoms.addBoxDerivatives( 3, Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1)) );
117 8888 : myatoms.addValue( 3, cp[1] );
118 :
119 8888 : addAtomDerivatives( 4, 0, crossProduct( Vector(0,0,-1.0), d2 ), myatoms );
120 8888 : if( myatoms.getNumberOfAtoms()==3 ) {
121 8888 : addAtomDerivatives( 4, 1, crossProduct( Vector(0,0,+1.0), d2 ) + crossProduct( Vector(0,0,-1.0), d1 ), myatoms );
122 8888 : addAtomDerivatives( 4, 2, crossProduct( Vector(0,0,+1.0), d1 ), myatoms);
123 : } else {
124 0 : addAtomDerivatives( 4, 1, crossProduct( Vector(0,0,-1.0), d2 ), myatoms);
125 0 : addAtomDerivatives( 4, 2, crossProduct( Vector(0,0,-1.0), d1 ), myatoms);
126 0 : addAtomDerivatives( 4, 3, crossProduct( Vector(0,0,+1.0), d1 ), myatoms);
127 : }
128 8888 : myatoms.addBoxDerivatives( 4, Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1)) );
129 8888 : myatoms.addValue( 4, cp[2] );
130 8888 : }
131 :
132 : // Vector MoleculePlane::getCentralAtom(){
133 : // Vector com; com.zero();
134 : // if( getNAtoms()==3 ){
135 : // com+=(1.0/3.0)*getPosition(0);
136 : // com+=(1.0/3.0)*getPosition(1);
137 : // com+=(1.0/3.0)*getPosition(2);
138 : // addCentralAtomDerivatives( 0, (1.0/3.0)*Tensor::identity() );
139 : // addCentralAtomDerivatives( 1, (1.0/3.0)*Tensor::identity() );
140 : // addCentralAtomDerivatives( 2, (1.0/3.0)*Tensor::identity() );
141 : // return com;
142 : // }
143 : // com+=0.25*getPosition(0);
144 : // com+=0.25*getPosition(1);
145 : // com+=0.25*getPosition(2);
146 : // com+=0.25*getPosition(3);
147 : // addCentralAtomDerivatives( 0, 0.25*Tensor::identity() );
148 : // addCentralAtomDerivatives( 1, 0.25*Tensor::identity() );
149 : // addCentralAtomDerivatives( 2, 0.25*Tensor::identity() );
150 : // addCentralAtomDerivatives( 3, 0.25*Tensor::identity() );
151 : // return com;
152 : // }
153 :
154 : }
155 : }
|
