Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2025 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 : #ifndef __PLUMED_secondarystructure_SecondaryStructureDRMSD_h
23 : #define __PLUMED_secondarystructure_SecondaryStructureDRMSD_h
24 : #include "SecondaryStructureBase.h"
25 : #include "tools/Matrix.h"
26 :
27 : #include<vector>
28 :
29 : namespace PLMD {
30 : namespace secondarystructure {
31 :
32 : ///This is a sort of compact std::vector<std::vector<T>>
33 : template<typename T>
34 : class flexibleMemory {
35 : ///Stores the data
36 : std::vector<T> data{};
37 : ///Stores the boundaries
38 : std::vector<size_t> sizes{};
39 : //helpers for openACC
40 : T* d;
41 : size_t* sz;
42 : void update() {
43 158 : d=data.data();
44 158 : sz=sizes.data();
45 66 : }
46 : public:
47 : flexibleMemory() {
48 : update();
49 : };
50 : ~flexibleMemory() =default;
51 : flexibleMemory(const flexibleMemory& m) :
52 : data(m.data),
53 : sizes(m.sizes) {
54 : update();
55 : }
56 66 : flexibleMemory& operator=(const flexibleMemory& m) {
57 66 : if (this!=&m) {
58 66 : data=m.data;
59 66 : sizes=m.sizes;
60 : update();
61 : }
62 66 : return *this;
63 : }
64 : constexpr size_t size() const {
65 : return sizes.size();
66 : }
67 : //Access the i-th view of the data, use this to access the elements
68 191232 : constexpr View< const T> get(size_t i) const {
69 : //use this at the start of your function, to not calculate again and again the initial index
70 : size_t init=0;
71 257236 : for(size_t j=0; j<i; j++) {
72 66004 : init+=sz[j];
73 : }
74 191232 : return View<const T>(d+init,sz[i]);
75 : }
76 92 : void extend(const std::vector<T>& v) {
77 92 : data.insert(data.end(),v.begin(),v.end());
78 92 : sizes.push_back(v.size());
79 : update();
80 92 : }
81 : void toACCDevice()const {
82 : #pragma acc enter data copyin(this[0:1], d[0:data.size()], sz[0:sizes.size()])
83 : }
84 : void removeFromACCDevice() const {
85 : #pragma acc exit data delete(sz[0:sizes.size()], d[0:data.size()], this[0:1])
86 : }
87 : };
88 :
89 : template <typename T>
90 : class SecondaryStructureDRMSDInput {
91 : private:
92 : struct pairDistance {
93 : T distance;
94 : unsigned i;
95 : unsigned j;
96 : };
97 : /// The list of reference configurations
98 : flexibleMemory<pairDistance> drmsd_targets{};
99 : flexibleMemory<unsigned> drmsd_atoms{};
100 : /// The general input for the secondary structure variable
101 : public:
102 148 : static void registerKeywords( Keywords& keys ) {
103 148 : keys.setDisplayName("SECONDARY_STRUCTURE_DRMSD");
104 148 : }
105 : static constexpr bool needsBondLength=true;
106 : typedef T precision;
107 : size_t natoms{0};
108 : size_t nstructures{0};
109 : size_t nindices_per_task{0};
110 : /// Are we operating without periodic boundary conditions
111 : bool nopbc{false};
112 : /// Variables for strands cutoff
113 : bool align_strands{false};
114 : /// The atoms involved in each of the secondary structure segments
115 : Matrix<unsigned> colvar_atoms{};
116 95616 : static void calculateDistance( const unsigned n,
117 : const bool noderiv,
118 : const SecondaryStructureDRMSDInput& actiondata,
119 : View<VectorT<precision>> pos,
120 : ColvarOutput<precision>& output ) {
121 95616 : const auto targetList = actiondata.drmsd_atoms.get(n);
122 : const auto targetAtoms=targetList.size();
123 95616 : if( !noderiv ) {
124 47244 : output.virial.set( n, TensorT<precision>(0,0,0,0,0,0,0,0,0) );
125 1463592 : for(unsigned i=0; i<targetAtoms; ++i ) {
126 1416348 : output.derivs[n][targetList[i]] =VectorT<precision>(0.0,0.0,0.0);
127 : }
128 : }
129 :
130 : double drmsd = 0.0;
131 95616 : Vector distance;
132 95616 : Vector dd;
133 :
134 95616 : const auto target = actiondata.drmsd_targets.get(n);
135 : const auto targetSize=target.size();
136 :
137 38981184 : for (unsigned i=0; i<targetSize; ++i) {
138 : const auto& [reference,k,j] = target[i];
139 38885568 : distance=delta( pos[k], pos[j] );
140 38885568 : const double len = distance.modulo();
141 38885568 : const double diff = len - reference;
142 38885568 : drmsd += diff*diff;
143 :
144 38885568 : if( !noderiv ) {
145 19213320 : dd=distance*(diff / len);
146 19213320 : output.derivs[n][k] -= dd;//-3 mul
147 19213320 : output.derivs[n][j] += dd;//-3mul
148 19213320 : output.virial.set( n, output.virial[n] - TensorT<precision>(dd,distance) );//-9 mul
149 : }
150 : }
151 95616 : const double inpairs = 1./static_cast<double>(targetSize);
152 95616 : output.values[n] = sqrt(inpairs*drmsd);
153 :
154 95616 : if( !noderiv ) {
155 : // The derivative of the square root is not defined when the structure is exactly the reference.
156 : // Using zero here avoids that 0*inf=nan is passed to the forces
157 47244 : const double scalef = output.values[n]>0.0 ? inpairs / output.values[n] : 0.0;
158 :
159 47244 : PLMD::LoopUnroller<9>::_mul(output.virial.getView(n).data(),scalef);
160 : // The atoms that do not take part in any pair are not in targetList: looping on
161 : // actiondata.natoms would read past the end of targetList
162 1463592 : for(unsigned i=0; i<targetAtoms; ++i ) {
163 1416348 : output.derivs[n][targetList[i]] *= scalef;
164 : }
165 : }
166 95616 : }
167 :
168 47 : void setReferenceStructure( const std::string& /*type*/, double bondlength, std::vector<Vector>& structure ) {
169 47 : std::vector<pairDistance> targets;
170 : //set is ordered and contains no duplicated data
171 : std::set<unsigned> atoms_targets;
172 : //targets.reserve( (structure.size()*(structure.size()-1))/2 );?
173 1330 : for(unsigned i=0; i+1<structure.size(); ++i) {
174 20435 : for(unsigned j=i+1; j<structure.size(); ++j) {
175 19152 : double distance = delta( structure[i], structure[j] ).modulo();
176 19152 : if(distance > bondlength) {
177 17906 : targets.emplace_back(pairDistance{distance,i,j});
178 17906 : atoms_targets.emplace(i);
179 17906 : atoms_targets.emplace(j);
180 : }
181 : }
182 : }
183 47 : if( targets.empty() ) {
184 2 : plumed_merror("all the distances in the reference structure are shorter than BONDLENGTH so there is nothing to compare");
185 : }
186 46 : drmsd_targets.extend( targets );
187 :
188 93 : drmsd_atoms.extend(std::vector<unsigned> {atoms_targets.begin(),atoms_targets.end()});
189 46 : nstructures = drmsd_targets.size();
190 46 : if (natoms==0) {
191 33 : natoms = structure.size();
192 13 : } else if (natoms!=structure.size()) {
193 0 : plumed_merror("input structures have a different number of atoms");
194 : }
195 :
196 46 : }
197 : void toACCDevice()const {
198 : #pragma acc enter data copyin(this[0:1], natoms,nstructures,nindices_per_task,nopbc,align_strands)
199 : drmsd_targets.toACCDevice();
200 : drmsd_atoms.toACCDevice();
201 : colvar_atoms.toACCDevice();
202 : }
203 : void removeFromACCDevice() const {
204 : colvar_atoms.removeFromACCDevice();
205 : drmsd_atoms.removeFromACCDevice();
206 : drmsd_targets.removeFromACCDevice();
207 : #pragma acc exit data delete(align_strands,nopbc,nindices_per_task,nstructures,natoms,this[0:1])
208 :
209 : }
210 : };
211 : }
212 : }
213 : #endif //__PLUMED_secondarystructure_SecondaryStructureDRMSD_h
|
