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 16 : void OrientationSphere::registerKeywords( Keywords& keys ) {
32 16 : multicolvar::MultiColvarBase::registerKeywords( keys );
33 32 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
34 32 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
35 32 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
36 32 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
37 32 : 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 48 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
42 48 : keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
43 64 : keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
44 32 : keys.use("LOWEST"); keys.use("HIGHEST");
45 16 : }
46 :
47 11 : OrientationSphere::OrientationSphere(const ActionOptions&ao):
48 : Action(ao),
49 11 : MultiColvarBase(ao)
50 : {
51 11 : 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 22 : std::string sw, errors; parse("SWITCH",sw);
54 11 : if(sw.length()>0) {
55 11 : 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 11 : log.printf(" degree of overlap in orientation between central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
64 22 : 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 11 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
67 11 : setLinkCellCutoff( switchingFunction.get_dmax() );
68 11 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms );
69 11 : }
70 :
71 7782 : double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
72 7782 : double sw, value=0, denom=0, dfunc; Vector ddistance;
73 7782 : unsigned ncomponents=getBaseMultiColvar(0)->getNumberOfQuantities();
74 7782 : std::vector<double> catom_orient( ncomponents ), this_orient( ncomponents );
75 7782 : std::vector<double> this_der( ncomponents ), catom_der( ncomponents );
76 :
77 7782 : getInputData( 0, true, myatoms, catom_orient );
78 7782 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
79 :
80 1379893 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
81 : Vector& distance=myatoms.getPosition(i);
82 : double d2;
83 2611747 : if ( (d2=distance[0]*distance[0])<rcut2 &&
84 1239636 : (d2+=distance[1]*distance[1])<rcut2 &&
85 2501268 : (d2+=distance[2]*distance[2])<rcut2 &&
86 : d2>epsilon ) {
87 :
88 1069658 : sw = switchingFunction.calculateSqr( d2, dfunc );
89 :
90 1069658 : getInputData( i, true, myatoms, this_orient );
91 : // Calculate the dot product wrt to this position
92 1069658 : double f_dot = computeVectorFunction( distance, catom_orient, this_orient, ddistance, catom_der, this_der );
93 :
94 1069658 : if( !doNotCalculateDerivatives() ) {
95 233625 : for(unsigned k=2; k<catom_orient.size(); ++k) { this_der[k]*=sw; catom_der[k]*=sw; }
96 10227 : MultiValue& myder1=getInputDerivatives( i, true, myatoms );
97 10227 : mergeInputDerivatives( 1, 2, this_orient.size(), 0, catom_der, myder0, myatoms );
98 10227 : mergeInputDerivatives( 1, 2, catom_der.size(), i, this_der, myder1, myatoms );
99 10227 : addAtomDerivatives( 1, 0, f_dot*(-dfunc)*distance - sw*ddistance, myatoms );
100 10227 : addAtomDerivatives( 1, i, f_dot*(dfunc)*distance + sw*ddistance, myatoms );
101 10227 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct(distance,ddistance) );
102 10227 : myder1.clearAll();
103 :
104 10227 : addAtomDerivatives( -1, 0, (-dfunc)*distance, myatoms );
105 10227 : addAtomDerivatives( -1, i, (dfunc)*distance, myatoms );
106 10227 : myatoms.addTemporyBoxDerivatives( (-dfunc)*Tensor(distance,distance) );
107 :
108 : }
109 1069658 : value += sw*f_dot;
110 1069658 : denom += sw;
111 : }
112 : }
113 7782 : double rdenom, df2, pref=calculateCoordinationPrefactor( denom, df2 );
114 7782 : if( std::fabs(denom)>epsilon ) { rdenom = 1.0 / denom; }
115 0 : else { plumed_assert(std::fabs(value)<epsilon); rdenom=1.0; }
116 :
117 : // Now divide everything
118 7782 : double rdenom2=rdenom*rdenom;
119 7782 : updateActiveAtoms( myatoms ); MultiValue& myvals=myatoms.getUnderlyingMultiValue();
120 57486 : for(unsigned i=0; i<myvals.getNumberActive(); ++i) {
121 49704 : unsigned ider=myvals.getActiveIndex(i);
122 : double dgd=myvals.getTemporyDerivative(ider);
123 49704 : myvals.setDerivative( 1, ider, rdenom*(pref*myvals.getDerivative(1,ider)+value*df2*dgd) - (value*pref*dgd)*rdenom2 );
124 : }
125 :
126 15564 : return pref*rdenom*value;
127 : }
128 :
129 : }
130 : }
131 :
|
