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 "OrientationSphere.h"
23 : #include "core/PlumedMain.h"
24 : #include "VectorMultiColvar.h"
25 :
26 : using namespace std;
27 :
28 : namespace PLMD {
29 : namespace crystallization {
30 :
31 31 : void OrientationSphere::registerKeywords( Keywords& keys ) {
32 31 : multicolvar::MultiColvarBase::registerKeywords( keys );
33 62 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
34 62 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
35 62 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
36 62 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
37 62 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
38 : "The following provides information on the \\ref switchingfunction that are available. "
39 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
40 : // Use actionWithDistributionKeywords
41 31 : keys.use("SPECIES");
42 31 : keys.use("SPECIESA");
43 31 : keys.use("SPECIESB");
44 31 : keys.use("MEAN");
45 31 : keys.use("MORE_THAN");
46 31 : keys.use("LESS_THAN");
47 31 : keys.use("MIN");
48 31 : keys.use("BETWEEN");
49 31 : keys.use("HISTOGRAM");
50 31 : keys.use("MOMENTS");
51 31 : keys.use("LOWEST");
52 31 : keys.use("HIGHEST");
53 31 : }
54 :
55 11 : OrientationSphere::OrientationSphere(const ActionOptions&ao):
56 : Action(ao),
57 11 : MultiColvarBase(ao) {
58 11 : if( getNumberOfBaseMultiColvars()>1 ) {
59 0 : warning("not sure if orientation sphere works with more than one base multicolvar - check numerical derivatives");
60 : }
61 : // Read in the switching function
62 : std::string sw, errors;
63 22 : parse("SWITCH",sw);
64 11 : if(sw.length()>0) {
65 11 : switchingFunction.set(sw,errors);
66 : } else {
67 0 : double r_0=-1.0, d_0;
68 : int nn, mm;
69 0 : parse("NN",nn);
70 0 : parse("MM",mm);
71 0 : parse("R_0",r_0);
72 0 : parse("D_0",d_0);
73 0 : if( r_0<0.0 ) {
74 0 : error("you must set a value for R_0");
75 : }
76 0 : switchingFunction.set(nn,mm,r_0,d_0);
77 : }
78 11 : log.printf(" degree of overlap in orientation between central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
79 22 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
80 : // Set the link cell cutoff
81 11 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
82 11 : setLinkCellCutoff( switchingFunction.get_dmax() );
83 : std::vector<AtomNumber> all_atoms;
84 11 : setupMultiColvarBase( all_atoms );
85 11 : }
86 :
87 7782 : double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
88 7782 : double sw, value=0, denom=0, dfunc;
89 7782 : Vector ddistance;
90 7782 : unsigned ncomponents=getBaseMultiColvar(0)->getNumberOfQuantities();
91 7782 : std::vector<double> catom_orient( ncomponents ), this_orient( ncomponents );
92 7782 : std::vector<double> this_der( ncomponents ), catom_der( ncomponents );
93 :
94 7782 : getInputData( 0, true, myatoms, catom_orient );
95 7782 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
96 :
97 1379893 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
98 : Vector& distance=myatoms.getPosition(i);
99 : double d2;
100 2611747 : if ( (d2=distance[0]*distance[0])<rcut2 &&
101 1239636 : (d2+=distance[1]*distance[1])<rcut2 &&
102 2501268 : (d2+=distance[2]*distance[2])<rcut2 &&
103 : d2>epsilon ) {
104 :
105 1069658 : sw = switchingFunction.calculateSqr( d2, dfunc );
106 :
107 1069658 : getInputData( i, true, myatoms, this_orient );
108 : // Calculate the dot product wrt to this position
109 1069658 : double f_dot = computeVectorFunction( distance, catom_orient, this_orient, ddistance, catom_der, this_der );
110 :
111 1069658 : if( !doNotCalculateDerivatives() ) {
112 233625 : for(unsigned k=2; k<catom_orient.size(); ++k) {
113 223398 : this_der[k]*=sw;
114 223398 : catom_der[k]*=sw;
115 : }
116 10227 : MultiValue& myder1=getInputDerivatives( i, true, myatoms );
117 10227 : mergeInputDerivatives( 1, 2, this_orient.size(), 0, catom_der, myder0, myatoms );
118 10227 : mergeInputDerivatives( 1, 2, catom_der.size(), i, this_der, myder1, myatoms );
119 10227 : addAtomDerivatives( 1, 0, f_dot*(-dfunc)*distance - sw*ddistance, myatoms );
120 10227 : addAtomDerivatives( 1, i, f_dot*(dfunc)*distance + sw*ddistance, myatoms );
121 10227 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct(distance,ddistance) );
122 10227 : myder1.clearAll();
123 :
124 10227 : addAtomDerivatives( -1, 0, (-dfunc)*distance, myatoms );
125 10227 : addAtomDerivatives( -1, i, (dfunc)*distance, myatoms );
126 10227 : myatoms.addTemporyBoxDerivatives( (-dfunc)*Tensor(distance,distance) );
127 :
128 : }
129 1069658 : value += sw*f_dot;
130 1069658 : denom += sw;
131 : }
132 : }
133 7782 : double rdenom, df2, pref=calculateCoordinationPrefactor( denom, df2 );
134 7782 : if( std::fabs(denom)>epsilon ) {
135 7782 : rdenom = 1.0 / denom;
136 : } else {
137 0 : plumed_assert(std::fabs(value)<epsilon);
138 : rdenom=1.0;
139 : }
140 :
141 : // Now divide everything
142 7782 : double rdenom2=rdenom*rdenom;
143 7782 : updateActiveAtoms( myatoms );
144 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
145 57486 : for(unsigned i=0; i<myvals.getNumberActive(); ++i) {
146 49704 : unsigned ider=myvals.getActiveIndex(i);
147 : double dgd=myvals.getTemporyDerivative(ider);
148 49704 : myvals.setDerivative( 1, ider, rdenom*(pref*myvals.getDerivative(1,ider)+value*df2*dgd) - (value*pref*dgd)*rdenom2 );
149 : }
150 :
151 15564 : return pref*rdenom*value;
152 : }
153 :
154 : }
155 : }
156 :
|
