Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-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 "tools/SwitchingFunction.h"
23 : #include "core/ActionRegister.h"
24 : #include "multicolvar/AtomValuePack.h"
25 : #include "VectorMultiColvar.h"
26 :
27 : //+PLUMEDOC MCOLVAR BOND_DIRECTIONS
28 : /*
29 : Calculate the vectors connecting atoms that are within cutoff defined using a switching function.
30 :
31 : \par Examples
32 :
33 : */
34 : //+ENDPLUMEDOC
35 :
36 : namespace PLMD {
37 : namespace crystallization {
38 :
39 : class BondOrientation : public VectorMultiColvar {
40 : private:
41 : double rcut2;
42 : SwitchingFunction switchingFunction;
43 : public:
44 : static void registerKeywords( Keywords& keys );
45 : explicit BondOrientation( const ActionOptions& ao );
46 : double calculateWeight( const unsigned& current, const double& weight, multicolvar::AtomValuePack& myatoms ) const override;
47 : void calculateVector( multicolvar::AtomValuePack& myatoms ) const override;
48 : };
49 :
50 13789 : PLUMED_REGISTER_ACTION(BondOrientation,"BOND_DIRECTIONS")
51 :
52 6 : void BondOrientation::registerKeywords( Keywords& keys ) {
53 6 : VectorMultiColvar::registerKeywords( keys );
54 12 : keys.add("numbered","ATOMS","the atoms involved in each of the vectors you wish to calculate. "
55 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one vector will be "
56 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
57 : "specify the indices of two atoms). The eventual number of quantities calculated by this "
58 : "action will depend on what functions of the distribution you choose to calculate.");
59 12 : keys.reset_style("ATOMS","atoms");
60 12 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
61 12 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
62 : "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
63 12 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
64 : "in GROUPB. This must be used in conjunction with GROUPA.");
65 12 : keys.add("compulsory","NN","12","The n parameter of the switching function ");
66 12 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
67 12 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
68 12 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
69 12 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
70 : "The following provides information on the \\ref switchingfunction that are available. "
71 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
72 6 : keys.use("VMEAN");
73 6 : keys.use("VSUM");
74 6 : }
75 :
76 2 : BondOrientation::BondOrientation( const ActionOptions& ao ):
77 : Action(ao),
78 2 : VectorMultiColvar(ao) {
79 : // Read atoms
80 2 : weightHasDerivatives=true;
81 : std::vector<AtomNumber> all_atoms;
82 4 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
83 2 : if( atom_lab.size()==0 ) {
84 0 : readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
85 : }
86 2 : setupMultiColvarBase( all_atoms );
87 : // Read in the switching function
88 : std::string sw, errors;
89 4 : parse("SWITCH",sw);
90 2 : if(sw.length()>0) {
91 2 : switchingFunction.set(sw,errors);
92 : } else {
93 0 : double r_0=-1.0, d_0;
94 : int nn, mm;
95 0 : parse("NN",nn);
96 0 : parse("MM",mm);
97 0 : parse("R_0",r_0);
98 0 : parse("D_0",d_0);
99 0 : if( r_0<0.0 ) {
100 0 : error("you must set a value for R_0");
101 : }
102 0 : switchingFunction.set(nn,mm,r_0,d_0);
103 : }
104 2 : log.printf(" orientation of those bonds with lengths are less than %s\n",( switchingFunction.description() ).c_str() );
105 : // Set the link cell cutoff
106 2 : setLinkCellCutoff( switchingFunction.get_dmax() );
107 2 : double rcut = switchingFunction.get_dmax();
108 2 : rcut2 = rcut*rcut;
109 : // Set the dimensionality of the vectors
110 2 : setVectorDimensionality(3);
111 2 : }
112 :
113 150 : double BondOrientation::calculateWeight( const unsigned& current, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
114 150 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
115 150 : double distm=distance.modulo2();
116 150 : if( distm>rcut2 ) {
117 : return 0.0;
118 : }
119 60 : double df, ww=switchingFunction.calculateSqr( distm, df );
120 : // Derivatives of weights
121 60 : addAtomDerivatives( 0, 0, -df*weight*distance, myatoms );
122 60 : addAtomDerivatives( 0, 1, df*weight*distance, myatoms );
123 60 : myatoms.addBoxDerivatives( 0, (-df)*weight*Tensor(distance,distance) );
124 60 : return ww;
125 : }
126 :
127 60 : void BondOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
128 60 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
129 :
130 60 : addAtomDerivatives( 2, 0, Vector(-1.0,0,0), myatoms );
131 60 : addAtomDerivatives( 2, 1, Vector(+1.0,0,0), myatoms );
132 60 : myatoms.addBoxDerivatives( 2, Tensor(distance,Vector(-1.0,0,0)) );
133 60 : myatoms.addValue( 2, distance[0] );
134 :
135 60 : addAtomDerivatives( 3, 0, Vector(0,-1.0,0), myatoms );
136 60 : addAtomDerivatives( 3, 1, Vector(0,+1.0,0), myatoms );
137 60 : myatoms.addBoxDerivatives( 3, Tensor(distance,Vector(0,-1.0,0)) );
138 60 : myatoms.addValue( 3, distance[1] );
139 :
140 60 : addAtomDerivatives( 4, 0, Vector(0,0,-1.0), myatoms );
141 60 : addAtomDerivatives( 4, 1, Vector(0,0,+1.0), myatoms );
142 60 : myatoms.addBoxDerivatives( 4, Tensor(distance,Vector(0,0,-1.0)) );
143 60 : myatoms.addValue( 4, distance[2] );
144 60 : }
145 :
146 : }
147 : }
|
