Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "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 14 : void OrientationSphere::registerKeywords( Keywords& keys ) {
32 14 : multicolvar::MultiColvarFunction::registerKeywords( keys );
33 14 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
34 14 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
35 14 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
36 14 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
37 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
38 : "The following provides information on the \\ref switchingfunction that are available. "
39 14 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
40 : // Use actionWithDistributionKeywords
41 14 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
42 14 : keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
43 14 : keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
44 14 : keys.use("LOWEST"); keys.use("HIGHEST"); keys.remove("DATA");
45 14 : }
46 :
47 10 : OrientationSphere::OrientationSphere(const ActionOptions&ao):
48 : Action(ao),
49 10 : MultiColvarFunction(ao)
50 : {
51 10 : if( getNumberOfBaseMultiColvars()>1 ) warning("not sure if orientation sphere works with more than one base multicolvar - check numerical derivatives");
52 : // Read in the switching function
53 20 : std::string sw, errors; parse("SWITCH",sw);
54 10 : if(sw.length()>0) {
55 10 : switchingFunction.set(sw,errors);
56 : } else {
57 0 : double r_0=-1.0, d_0; int nn, mm;
58 0 : parse("NN",nn); parse("MM",mm);
59 0 : parse("R_0",r_0); parse("D_0",d_0);
60 0 : if( r_0<0.0 ) error("you must set a value for R_0");
61 0 : switchingFunction.set(nn,mm,r_0,d_0);
62 : }
63 10 : log.printf(" degree of overlap in orientation between central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
64 10 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
65 : // Set the link cell cutoff
66 10 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
67 10 : setLinkCellCutoff( switchingFunction.get_dmax() );
68 20 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms );
69 10 : }
70 :
71 7694 : double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
72 7694 : double d2, sw, value=0, denom=0, dfunc; Vector ddistance;
73 7694 : unsigned ncomponents=getBaseMultiColvar(0)->getNumberOfQuantities();
74 15388 : std::vector<double> catom_orient( ncomponents ), this_orient( ncomponents );
75 15388 : std::vector<double> this_der( ncomponents ), catom_der( ncomponents );
76 :
77 7694 : getInputData( 0, true, myatoms, catom_orient );
78 15388 : multicolvar::CatomPack atom0;
79 7694 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
80 7694 : if( !doNotCalculateDerivatives() ) atom0=getCentralAtomPackFromInput( myatoms.getIndex(0) );
81 :
82 1379189 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
83 1371495 : Vector& distance=myatoms.getPosition(i);
84 6592525 : if ( (d2=distance[0]*distance[0])<rcut2 &&
85 4735122 : (d2+=distance[1]*distance[1])<rcut2 &&
86 6069207 : (d2+=distance[2]*distance[2])<rcut2 &&
87 1069135 : d2>epsilon ) {
88 :
89 1069042 : sw = switchingFunction.calculateSqr( d2, dfunc );
90 :
91 1069042 : getInputData( i, true, myatoms, this_orient );
92 : // Calculate the dot product wrt to this position
93 1069042 : double f_dot = computeVectorFunction( distance, catom_orient, this_orient, ddistance, catom_der, this_der );
94 :
95 1069042 : if( !doNotCalculateDerivatives() ) {
96 9611 : for(unsigned k=2; k<catom_orient.size(); ++k) { this_der[k]*=sw; catom_der[k]*=sw; }
97 9611 : MultiValue& myder1=getInputDerivatives( i, true, myatoms );
98 9611 : mergeInputDerivatives( 1, 2, this_orient.size(), 0, catom_der, myder0, myatoms );
99 9611 : mergeInputDerivatives( 1, 2, catom_der.size(), i, this_der, myder1, myatoms );
100 9611 : myatoms.addComDerivatives( 1, f_dot*(-dfunc)*distance - sw*ddistance, atom0 );
101 9611 : addAtomDerivatives( 1, i, f_dot*(dfunc)*distance + sw*ddistance, myatoms );
102 9611 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct(distance,ddistance) );
103 9611 : myder1.clearAll();
104 :
105 9611 : myatoms.addComDerivatives( -1, (-dfunc)*distance, atom0 );
106 9611 : addAtomDerivatives( -1, i, (dfunc)*distance, myatoms );
107 9611 : myatoms.addTemporyBoxDerivatives( (-dfunc)*Tensor(distance,distance) );
108 :
109 : }
110 1069042 : value += sw*f_dot;
111 1069042 : denom += sw;
112 : }
113 : }
114 7694 : double rdenom, df2, pref=calculateCoordinationPrefactor( denom, df2 );
115 7694 : if( fabs(denom)>epsilon ) { rdenom = 1.0 / denom; }
116 0 : else { plumed_assert(fabs(value)<epsilon); rdenom=1.0; }
117 :
118 : // Now divide everything
119 7694 : double rdenom2=rdenom*rdenom;
120 7694 : updateActiveAtoms( myatoms ); MultiValue& myvals=myatoms.getUnderlyingMultiValue();
121 50270 : for(unsigned i=0; i<myvals.getNumberActive(); ++i) {
122 42576 : unsigned ider=myvals.getActiveIndex(i);
123 42576 : double dgd=myvals.getTemporyDerivative(ider);
124 42576 : myvals.setDerivative( 1, ider, rdenom*(pref*myvals.getDerivative(1,ider)+value*df2*dgd) - (value*pref*dgd)*rdenom2 );
125 : }
126 :
127 15388 : return pref*rdenom*value;
128 : }
129 :
130 : }
131 2523 : }
132 :
|
