Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "CubicHarmonicBase.h"
23 : #include "tools/SwitchingFunction.h"
24 :
25 : #include <string>
26 : #include <cmath>
27 :
28 : using namespace std;
29 :
30 : namespace PLMD {
31 : namespace crystallization {
32 :
33 18 : void CubicHarmonicBase::registerKeywords( Keywords& keys ) {
34 18 : multicolvar::MultiColvarBase::registerKeywords( keys );
35 18 : keys.use("SPECIES");
36 18 : keys.use("SPECIESA");
37 18 : keys.use("SPECIESB");
38 36 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
39 36 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
40 36 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
41 36 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
42 36 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
43 : "The following provides information on the \\ref switchingfunction that are available. "
44 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
45 36 : keys.add("compulsory","PHI","0.0","The Euler rotational angle phi");
46 36 : keys.add("compulsory","THETA","0.0","The Euler rotational angle theta");
47 36 : keys.add("compulsory","PSI","0.0","The Euler rotational angle psi");
48 36 : keys.addFlag("UNORMALIZED",false,"calculate the sum of the components of the vector rather than the mean");
49 : // Use actionWithDistributionKeywords
50 18 : keys.use("MEAN");
51 18 : keys.use("MORE_THAN");
52 18 : keys.use("LESS_THAN");
53 18 : keys.use("MAX");
54 18 : keys.use("MIN");
55 18 : keys.use("BETWEEN");
56 18 : keys.use("HISTOGRAM");
57 18 : keys.use("MOMENTS");
58 18 : keys.use("ALT_MIN");
59 18 : keys.use("LOWEST");
60 18 : keys.use("HIGHEST");
61 18 : }
62 :
63 6 : CubicHarmonicBase::CubicHarmonicBase(const ActionOptions&ao):
64 : Action(ao),
65 6 : MultiColvarBase(ao) {
66 : // Read in the switching function
67 : std::string sw, errors;
68 12 : parse("SWITCH",sw);
69 6 : if(sw.length()>0) {
70 6 : switchingFunction.set(sw,errors);
71 6 : if( errors.length()!=0 ) {
72 0 : error("problem reading SWITCH keyword : " + errors );
73 : }
74 : } else {
75 0 : double r_0=-1.0, d_0;
76 : int nn, mm;
77 0 : parse("NN",nn);
78 0 : parse("MM",mm);
79 0 : parse("R_0",r_0);
80 0 : parse("D_0",d_0);
81 0 : if( r_0<0.0 ) {
82 0 : error("you must set a value for R_0");
83 : }
84 0 : switchingFunction.set(nn,mm,r_0,d_0);
85 : }
86 :
87 : double phi, theta, psi;
88 6 : parse("PHI",phi);
89 6 : parse("THETA",theta);
90 6 : parse("PSI",psi);
91 6 : log.printf(" creating rotation matrix with Euler angles phi=%f, theta=%f and psi=%f\n",phi,theta,psi);
92 : // Calculate the rotation matrix http://mathworld.wolfram.com/EulerAngles.html
93 6 : rotationmatrix[0][0]=std::cos(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::sin(psi);
94 6 : rotationmatrix[0][1]=std::cos(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::sin(psi);
95 6 : rotationmatrix[0][2]=std::sin(psi)*std::sin(theta);
96 :
97 6 : rotationmatrix[1][0]=-std::sin(psi)*std::cos(phi)-std::cos(theta)*std::sin(phi)*std::cos(psi);
98 6 : rotationmatrix[1][1]=-std::sin(psi)*std::sin(phi)+std::cos(theta)*std::cos(phi)*std::cos(psi);
99 6 : rotationmatrix[1][2]=std::cos(psi)*std::sin(theta);
100 :
101 6 : rotationmatrix[2][0]=std::sin(theta)*std::sin(phi);
102 6 : rotationmatrix[2][1]=-std::sin(theta)*std::cos(phi);
103 6 : rotationmatrix[2][2]=std::cos(theta);
104 :
105 :
106 6 : log.printf(" measure crystallinity around central atom. Includes those atoms within %s\n",( switchingFunction.description() ).c_str() );
107 6 : parseFlag("UNORMALIZED",unormalized);
108 6 : if( unormalized ) {
109 0 : log.printf(" output sum of vector functions \n");
110 : } else {
111 6 : log.printf(" output mean of vector functions \n");
112 : }
113 : // Set the link cell cutoff
114 6 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
115 6 : setLinkCellCutoff( switchingFunction.get_dmax() );
116 : // And setup the ActionWithVessel
117 : std::vector<AtomNumber> all_atoms;
118 6 : setupMultiColvarBase( all_atoms );
119 6 : }
120 :
121 26142 : double CubicHarmonicBase::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
122 : double dfunc;
123 26142 : Vector rotatedis;
124 :
125 : // Calculate the coordination number
126 26142 : Vector myder, rotateder, fder;
127 : unsigned nat=myatoms.getNumberOfAtoms();
128 :
129 2517483 : for(unsigned i=1; i<nat; ++i) {
130 : Vector& distance=myatoms.getPosition(i);
131 :
132 : double d2;
133 3953726 : if ( (d2=distance[0]*distance[0])<rcut2 &&
134 1462385 : (d2+=distance[1]*distance[1])<rcut2 &&
135 3193615 : (d2+=distance[2]*distance[2])<rcut2 &&
136 : d2>epsilon ) {
137 :
138 328678 : double sw = switchingFunction.calculateSqr( d2, dfunc );
139 :
140 328678 : rotatedis[0]=rotationmatrix[0][0]*distance[0]
141 328678 : +rotationmatrix[0][1]*distance[1]
142 328678 : +rotationmatrix[0][2]*distance[2];
143 328678 : rotatedis[1]=rotationmatrix[1][0]*distance[0]
144 328678 : +rotationmatrix[1][1]*distance[1]
145 328678 : +rotationmatrix[1][2]*distance[2];
146 328678 : rotatedis[2]=rotationmatrix[2][0]*distance[0]
147 328678 : +rotationmatrix[2][1]*distance[1]
148 328678 : +rotationmatrix[2][2]*distance[2];
149 :
150 328678 : double tmp = calculateCubicHarmonic( rotatedis, d2, rotateder );
151 :
152 328678 : myder[0]=rotationmatrix[0][0]*rotateder[0]
153 328678 : +rotationmatrix[1][0]*rotateder[1]
154 328678 : +rotationmatrix[2][0]*rotateder[2];
155 328678 : myder[1]=rotationmatrix[0][1]*rotateder[0]
156 328678 : +rotationmatrix[1][1]*rotateder[1]
157 328678 : +rotationmatrix[2][1]*rotateder[2];
158 328678 : myder[2]=rotationmatrix[0][2]*rotateder[0]
159 328678 : +rotationmatrix[1][2]*rotateder[1]
160 328678 : +rotationmatrix[2][2]*rotateder[2];
161 :
162 328678 : fder = (+dfunc)*tmp*distance + sw*myder;
163 :
164 328678 : accumulateSymmetryFunction( 1, i, sw*tmp, fder, Tensor(distance,-fder), myatoms );
165 328678 : accumulateSymmetryFunction( -1, i, sw, (+dfunc)*distance, (-dfunc)*Tensor(distance,distance), myatoms );
166 : }
167 : }
168 : // values -> der of... value [0], weight[1], x coord [2], y, z... [more magic]
169 26142 : updateActiveAtoms( myatoms );
170 26142 : if( !unormalized ) {
171 26142 : myatoms.getUnderlyingMultiValue().quotientRule( 1, 1 );
172 : }
173 26142 : return myatoms.getValue(1); // this is equivalent to getting an "atomic" CV
174 : }
175 :
176 : }
177 : }
178 :
|
