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