Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2022.
3 :
4 : CVs originally developed by the Jochen Hub group from the University of Saarland (Germany)
5 : and adapted and implemented in PLUMED by Ary Lautaro Di Bartolo and Diego Masone from the
6 : National University of Cuyo (Argentina).
7 :
8 : The membranefusion module is free software: you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation, either version 3 of the License, or
11 : (at your option) any later version.
12 :
13 : The membranefusion module is distributed in the hope that it will be useful,
14 : but WITHOUT ANY WARRANTY; without even the implied warranty of
15 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 : GNU Lesser General Public License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
20 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
21 : #include "colvar/Colvar.h"
22 : #include "core/ActionRegister.h"
23 : #include <cmath>
24 : #ifdef _OPENMP
25 : #if _OPENMP >= 201307
26 : #include <omp.h>
27 : #endif
28 : #endif
29 :
30 : namespace PLMD {
31 : namespace membranefusion {
32 : //+PLUMEDOC MEMBRANEFUSIONMOD_COLVAR MEMFUSIONP
33 : /*
34 : Calculate a CV that can induce the formation of the hemifusion stalk between two initially flat and planar bilayers.
35 :
36 : Calculate the collective variable designed by Hub and collaborators \cite Hub2017 and
37 : implemented into PLUMED by Masone and collaborators \cite DiBartolo2022 .
38 : This CV is capable of inducing the formation of the hemifusion stalk between two initially flat and planar bilayers
39 : surrounded by water molecules.
40 :
41 : \f[
42 : \xi_f = \frac{1}{N_{sf}} \sum_{s=0}^{N_{sf}-1} \delta_{sf} (N_{sf}^{(p)})
43 : \f]
44 :
45 : Where \f$\xi_f\f$ is the CV, \f$N_{sf}\f$ is the number of slices of the cylinder that make up the CV,
46 : \f$\delta_{sf}\f$ is a continuos function in the interval [0 1] (\f$\delta_{sf} = 0\f$ for no beads in the slice s, and
47 : \f$\delta_{sf} = 1\f$ for 1 or more beads in the slice s) and \f$N_{sf}^{(p)}\f$ accounts for the number of tail beads
48 : within the slice s.
49 :
50 : \par Examples
51 :
52 : This example induces a hemifusion stalk (\f$\xi_f = 0.85\f$) from a pair of initially flat membranes (\f$\xi_f = 0.2\f$).
53 :
54 : \plumedfile
55 : lMem: GROUP ATOMS=1-12288 #All the lower membrane beads.
56 : uMem: GROUP ATOMS=12289-24576 #All the upper membrane beads.
57 : tails: GROUP ATOMS=8-24572:12,12-24576:12 #All the lipid tails beads (from the lower and upper membrane).
58 :
59 : memFusion: MEMFUSIONP UMEMBRANE=uMem LMEMBRANE=lMem TAILS=tails NSMEM=70 DSMEM=0.1 HMEM=0.25 RCYLMEM=1.75 ZETAMEM=0.5
60 :
61 : MOVINGRESTRAINT ...
62 : ARG=memFusion
63 : STEP0=0 AT0=0.2 KAPPA0=10000.0
64 : STEP1=500000 AT1=0.85 KAPPA1=10000.0
65 : ...
66 :
67 : PRINT ARG=memFusion FILE=COLVAR STRIDE=1
68 :
69 : \endplumedfile
70 :
71 : You can test this CV with another example in this <a href="https://github.com/lautarodibartolo/MemFusion/tree/main/ExampleParallel">GitHub folder</a>.
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 : class memFusionP : public Colvar {
77 : std::vector<AtomNumber> UMEM, LMEM, TAILS;
78 : std::vector<double> NSMEM, DSMEM, HMEM, RCYLMEM, ZETAMEM, ONEOVERS2C2CUTOFF, XCYL, YCYL;
79 :
80 : public:
81 : explicit memFusionP(const ActionOptions &);
82 : void calculate() override;
83 : static void registerKeywords(Keywords &keys);
84 : };
85 :
86 13787 : PLUMED_REGISTER_ACTION(memFusionP, "MEMFUSIONP")
87 :
88 5 : void memFusionP::registerKeywords(Keywords &keys) {
89 5 : Colvar::registerKeywords(keys);
90 10 : keys.add("atoms", "UMEMBRANE", "all the beads of the upper membrane");
91 10 : keys.add("atoms", "LMEMBRANE", "all the beads of the lower membrane");
92 10 : keys.add("atoms", "TAILS", "all the tail beads of the system");
93 10 : keys.add("compulsory", "NSMEM", "the number of slices of the membrane fusion cylinder in such a way that when the bilayers are flat and parallel the CV is equal to 0.2.");
94 10 : keys.add("optional", "DSMEM", "( default=0.1) thickness of the slices of the membrane fusion cylinder.");
95 10 : keys.add("optional", "HMEM", "( default=0.25 ) parameter of the step function θ(x,h) for the membrane fusion.");
96 10 : keys.add("optional", "RCYLMEM", "( default=1.75 ) the radius of the membrane fusion cylinder.");
97 10 : keys.add("optional", "ZETAMEM", "( default=0.5 ) occupation factor.");
98 10 : keys.add("optional", "ONEOVERS2C2CUTOFF", "( default=500 ) cut off large values for the derivative of the atan2 function.");
99 10 : keys.add("optional", "XCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
100 10 : keys.add("optional", "YCYL", "Y coordinate of the fixed cylinder, if not present this will be calculated.");
101 5 : }
102 :
103 1 : memFusionP::memFusionP(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao) {
104 2 : parseAtomList("UMEMBRANE", UMEM);
105 1 : if (UMEM.size() == 0) {
106 0 : error("UMEMBRANE has not any atom specified.");
107 : }
108 :
109 2 : parseAtomList("LMEMBRANE", LMEM);
110 1 : if (LMEM.size() == 0) {
111 0 : error("LMEMBRANE has not any atom specified.");
112 : }
113 :
114 2 : parseAtomList("TAILS", TAILS);
115 1 : if (TAILS.size() == 0) {
116 0 : error("TAILS has not any atom specified.");
117 : }
118 :
119 2 : parseVector("NSMEM", NSMEM);
120 1 : if (NSMEM.size() > 1) {
121 0 : error("NSMEM cannot take more than one value.");
122 : }
123 :
124 2 : parseVector("DSMEM", DSMEM);
125 1 : if (DSMEM.size() > 1) {
126 0 : error("DSMEM cannot take more than one value.");
127 : }
128 1 : if (DSMEM.size() == 0) {
129 0 : DSMEM.push_back(0.1);
130 : }
131 :
132 2 : parseVector("HMEM", HMEM);
133 1 : if (HMEM.size() > 1) {
134 0 : error("HMEM cannot take more than one value.");
135 : }
136 1 : if (HMEM.size() == 0) {
137 0 : HMEM.push_back(0.25);
138 : }
139 :
140 2 : parseVector("RCYLMEM", RCYLMEM);
141 1 : if (RCYLMEM.size() > 1) {
142 0 : error("RCYLMEM cannot take more than one value.");
143 : }
144 1 : if (RCYLMEM.size() == 0) {
145 0 : RCYLMEM.push_back(1.75);
146 : }
147 :
148 2 : parseVector("ZETAMEM", ZETAMEM);
149 1 : if (ZETAMEM.size() > 1) {
150 0 : error("ZETA cannot take more than one value.");
151 : }
152 1 : if (ZETAMEM.size() == 0) {
153 0 : ZETAMEM.push_back(0.5);
154 : }
155 :
156 2 : parseVector("ONEOVERS2C2CUTOFF", ONEOVERS2C2CUTOFF);
157 1 : if (ONEOVERS2C2CUTOFF.size() > 1) {
158 0 : error("ONEOVERS2C2CUTOFF cannot take more than one value.");
159 : }
160 1 : if (ONEOVERS2C2CUTOFF.size() == 0) {
161 1 : ONEOVERS2C2CUTOFF.push_back(500);
162 : }
163 :
164 2 : parseVector("XCYL", XCYL);
165 1 : if (XCYL.size() > 1) {
166 0 : error("XCYL cannot take more than one value.");
167 : }
168 1 : if (XCYL.size() == 0) {
169 1 : XCYL.push_back(-1.0);
170 : }
171 :
172 2 : parseVector("YCYL", YCYL);
173 1 : if (YCYL.size() > 1) {
174 0 : error("YCYL cannot take more than one value.");
175 : }
176 1 : if (YCYL.size() == 0) {
177 1 : YCYL.push_back(-1.0);
178 : }
179 :
180 1 : checkRead();
181 :
182 : std::vector<AtomNumber> atoms;
183 12289 : for (unsigned i = 0; i < UMEM.size(); i++) {
184 12288 : atoms.push_back(UMEM[i]);
185 : }
186 12289 : for (unsigned i = 0; i < LMEM.size(); i++) {
187 12288 : atoms.push_back(LMEM[i]);
188 : }
189 4097 : for (unsigned i = 0; i < TAILS.size(); i++) {
190 4096 : atoms.push_back(TAILS[i]);
191 : }
192 :
193 1 : addValueWithDerivatives();
194 1 : setNotPeriodic();
195 1 : requestAtoms(atoms);
196 1 : }
197 :
198 3 : void memFusionP::calculate() {
199 : /**************************
200 : * *
201 : * System *
202 : * *
203 : **************************/
204 :
205 : // Box dimensions.
206 3 : double Lx = getBox()[0][0], Ly = getBox()[1][1], Lz = getBox()[2][2];
207 :
208 : // Z center of the upper membrane (uMem) and lower membrane (lMem) for systems with PBC: https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions .
209 : double ZuMem, ZuMemcos = 0.0, ZuMemsin = 0.0, uMemAngle, ZlMem, ZlMemcos = 0.0, ZlMemsin = 0.0, lMemAngle;
210 :
211 : #ifdef _OPENMP
212 : #if _OPENMP >= 201307
213 3 : #pragma omp parallel for private(uMemAngle, lMemAngle) reduction(+:ZuMemcos, ZuMemsin, ZlMemcos, ZlMemsin)
214 : #endif
215 : #endif
216 : for (unsigned i = 0; i < UMEM.size(); i++) {
217 : uMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
218 : lMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + UMEM.size())))[2];
219 : ZuMemcos += cos(uMemAngle);
220 : ZuMemsin += sin(uMemAngle);
221 : ZlMemcos += cos(lMemAngle);
222 : ZlMemsin += sin(lMemAngle);
223 : }
224 :
225 3 : ZuMemcos = ZuMemcos / UMEM.size();
226 3 : ZuMemsin = ZuMemsin / UMEM.size();
227 3 : ZlMemcos = ZlMemcos / UMEM.size();
228 3 : ZlMemsin = ZlMemsin / UMEM.size();
229 :
230 3 : ZuMem = Lz * (atan2(-ZuMemsin, -ZuMemcos) + M_PI) / (2.0 * M_PI);
231 3 : ZlMem = Lz * (atan2(-ZlMemsin, -ZlMemcos) + M_PI) / (2.0 * M_PI);
232 :
233 : // Z center of the boths membranes (upper and lower).
234 3 : double ZMems = (ZuMem + ZlMem) / 2.0;
235 :
236 : /*************************
237 : * *
238 : * Xi_Mem *
239 : * *
240 : **************************/
241 :
242 : // Quantity of beads of the membranes.
243 3 : unsigned membraneBeads = UMEM.size() + LMEM.size();
244 :
245 : // Z distance from the lipid tail to the geometric center of both membranes.
246 : double ZTailDistance;
247 :
248 : // Z position of the first slice.
249 3 : double firstSliceZ_Mem = ZMems + (0.0 + 0.5 - NSMEM[0] / 2.0) * DSMEM[0];
250 :
251 : // Z distance between the first slice and the Z center of the membrane.
252 3 : double firstSliceZDist_Mem = pbcDistance(Vector(0.0, 0.0, firstSliceZ_Mem), Vector(0.0, 0.0, ZMems))[2];
253 :
254 : // Position in the cylinder.
255 : double PositionS_Mem;
256 :
257 : // Slices to analyze per particle.
258 3 : std::vector<unsigned> s1_Mem(TAILS.size()), s2_Mem(TAILS.size());
259 :
260 : // Mark the particles to analyze.
261 3 : std::vector<double> analyzeThisParticle_Mem(TAILS.size());
262 :
263 : // Eq. 7 Hub & Awasthi JCTC 2017.
264 3 : std::vector<double> faxial_Mem(TAILS.size() * NSMEM[0]);
265 :
266 : // Eq. 16 Hub & Awasthi JCTC 2017.
267 3 : std::vector<double> d_faxial_Mem_dz(TAILS.size() * NSMEM[0]);
268 :
269 : // Eq. 10 Hub & Awasthi JCTC 2017.
270 3 : std::vector<double> Fs_Mem(NSMEM[0]);
271 :
272 : // Eq. 11 Hub & Awasthi JCTC 2017.
273 3 : std::vector<double> ws_Mem(NSMEM[0]);
274 :
275 : // Eq. 10 Hub & Awasthi JCTC 2017.
276 : double W_Mem = 0.0;
277 :
278 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
279 3 : std::vector<double> sx_Mem(NSMEM[0]), sy_Mem(NSMEM[0]), cx_Mem(NSMEM[0]), cy_Mem(NSMEM[0]);
280 :
281 : // Eq. 10 Hub & Awasthi JCTC 2017.
282 : double Xsc_Mem = 0.0, Xcc_Mem = 0.0, Ysc_Mem = 0.0, Ycc_Mem = 0.0;
283 :
284 : // Aux.
285 : double x, aux;
286 :
287 : // Scaled position of the lipid tail respect the origin of coordinates.
288 3 : Vector TailPosition;
289 :
290 : // Thanks stack overflow.
291 : #ifdef _OPENMP
292 : #if _OPENMP >= 201307
293 : #pragma omp declare reduction(vec_double_plus : std::vector<double> : \
294 : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
295 : initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
296 : #endif
297 : #endif
298 :
299 : #ifdef _OPENMP
300 : #if _OPENMP >= 201307
301 3 : #pragma omp parallel for private(ZTailDistance, PositionS_Mem, TailPosition, x, aux) reduction(vec_double_plus:Fs_Mem, sx_Mem, sy_Mem, cx_Mem, cy_Mem)
302 : #endif
303 : #endif
304 : for (unsigned i = 0; i < TAILS.size(); i++) {
305 : ZTailDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + membraneBeads))[2];
306 : PositionS_Mem = (ZTailDistance + firstSliceZDist_Mem) / DSMEM[0];
307 : // If the following condition is met the particle is in the Z space of the cylinder.
308 : if ((PositionS_Mem >= (-0.5 - HMEM[0])) && (PositionS_Mem <= (NSMEM[0] + 0.5 - 1.0 + HMEM[0]))) {
309 : analyzeThisParticle_Mem[i] = 1.0;
310 : // Defining the slices to analyze each particle.
311 : if (PositionS_Mem < 1) {
312 : s1_Mem[i] = 0;
313 : s2_Mem[i] = 2;
314 : } else if (PositionS_Mem <= (NSMEM[0] - 2.0)) {
315 : s1_Mem[i] = floor(PositionS_Mem) - 1;
316 : s2_Mem[i] = floor(PositionS_Mem) + 1;
317 : } else {
318 : s1_Mem[i] = NSMEM[0] - 3;
319 : s2_Mem[i] = NSMEM[0] - 1;
320 : }
321 :
322 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
323 :
324 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
325 : x = (ZTailDistance - (s + 0.5 - NSMEM[0] / 2.0) * DSMEM[0]) * 2.0 / DSMEM[0];
326 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0]))) {
327 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0]))) {
328 : faxial_Mem[i + TAILS.size() * s] = 1.0;
329 : Fs_Mem[s] += 1.0;
330 : sx_Mem[s] += sin(2.0 * M_PI * TailPosition[0]);
331 : sy_Mem[s] += sin(2.0 * M_PI * TailPosition[1]);
332 : cx_Mem[s] += cos(2.0 * M_PI * TailPosition[0]);
333 : cy_Mem[s] += cos(2.0 * M_PI * TailPosition[1]);
334 : } else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0]))) {
335 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
336 : faxial_Mem[i + TAILS.size() * s] = aux;
337 : d_faxial_Mem_dz[i + TAILS.size() * s] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * 2.0 / DSMEM[0];
338 : Fs_Mem[s] += aux;
339 : sx_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[0]);
340 : sy_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[1]);
341 : cx_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[0]);
342 : cy_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[1]);
343 : } else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0]))) {
344 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
345 : faxial_Mem[i + TAILS.size() * s] = aux;
346 : d_faxial_Mem_dz[i + TAILS.size() * s] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * 2.0 / DSMEM[0];
347 : Fs_Mem[s] += aux;
348 : sx_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[0]));
349 : sy_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[1]));
350 : cx_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[0]));
351 : cy_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[1]));
352 : }
353 : }
354 : }
355 : }
356 : }
357 :
358 213 : for (unsigned s = 0; s < NSMEM[0]; s++) {
359 210 : if (Fs_Mem[s] != 0.0) {
360 106 : ws_Mem[s] = tanh(Fs_Mem[s]);
361 106 : W_Mem += ws_Mem[s];
362 106 : sx_Mem[s] = sx_Mem[s] / Fs_Mem[s];
363 106 : sy_Mem[s] = sy_Mem[s] / Fs_Mem[s];
364 106 : cx_Mem[s] = cx_Mem[s] / Fs_Mem[s];
365 106 : cy_Mem[s] = cy_Mem[s] / Fs_Mem[s];
366 106 : Xsc_Mem += sx_Mem[s] * ws_Mem[s];
367 106 : Ysc_Mem += sy_Mem[s] * ws_Mem[s];
368 106 : Xcc_Mem += cx_Mem[s] * ws_Mem[s];
369 106 : Ycc_Mem += cy_Mem[s] * ws_Mem[s];
370 : }
371 : }
372 :
373 3 : Xsc_Mem = Xsc_Mem / W_Mem;
374 3 : Ysc_Mem = Ysc_Mem / W_Mem;
375 3 : Xcc_Mem = Xcc_Mem / W_Mem;
376 3 : Ycc_Mem = Ycc_Mem / W_Mem;
377 :
378 : // Eq. 12 Hub & Awasthi JCTC 2017.
379 : double Xcyl_Mem, Ycyl_Mem;
380 :
381 3 : if ((XCYL[0] > 0.0) && (YCYL[0] > 0.0)) {
382 : Xcyl_Mem = XCYL[0];
383 : Ycyl_Mem = YCYL[0];
384 : } else {
385 3 : Xcyl_Mem = (atan2(-Xsc_Mem, -Xcc_Mem) + M_PI) * Lx / (2 * M_PI);
386 3 : Ycyl_Mem = (atan2(-Ysc_Mem, -Ycc_Mem) + M_PI) * Ly / (2 * M_PI);
387 : }
388 :
389 : // Eq. 25, 26 and 27 Hub & Awasthi JCTC 2017.
390 : double d_sx_Mem_dx, d_sx_Mem_dz, d_sy_Mem_dy, d_sy_Mem_dz, d_cx_Mem_dx, d_cx_Mem_dz, d_cy_Mem_dy, d_cy_Mem_dz;
391 :
392 : // Eq. 29 Hub & Awasthi JCTC 2017.
393 : double d_ws_Mem_dz;
394 :
395 : // Eq. 31, 32 and 33 Hub & Awasthi JCTC 2017
396 : double d_Xsc_Mem_dx, d_Xsc_Mem_dz, d_Xcc_Mem_dx, d_Xcc_Mem_dz, d_Ysc_Mem_dy, d_Ysc_Mem_dz, d_Ycc_Mem_dy, d_Ycc_Mem_dz;
397 :
398 : // Center of the cylinder. XY components are calculated (or defined), Z is the Z geometric center of the membranes of the system.
399 3 : Vector xyzCyl_Mem = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMems));
400 :
401 : // Distances from the lipid tails to center of the cylinder.
402 3 : std::vector<Vector> CylDistances_Mem(TAILS.size());
403 :
404 : // XY distance from the lipid tails to the center of the cylinder.
405 : double ri_Mem;
406 :
407 : // Eq. 8 Hub & Awasthi JCTC 2017.
408 : double fradial_Mem = 0;
409 :
410 : // Eq. 15 Hub & Awasthi JCTC 2017.
411 3 : std::vector<double> d_fradial_Mem_dx(TAILS.size()), d_fradial_Mem_dy(TAILS.size());
412 :
413 : // Eq. 35, 36, 37 and 38 Hub & Awasthi JCTC 2017.
414 3 : std::vector<double> d_Xcyl_Mem_dx(TAILS.size()), d_Xcyl_Mem_dz(TAILS.size()), d_Ycyl_Mem_dy(TAILS.size()), d_Ycyl_Mem_dz(TAILS.size());
415 :
416 : // To avoid rare instabilities auxX_Mem and auxY_Mem are truncated at a configurable value (default = 500).
417 3 : double auxX_Mem = (1 / (pow(Xsc_Mem, 2) + pow(Xcc_Mem, 2))), auxY_Mem = (1 / (pow(Ysc_Mem, 2) + pow(Ycc_Mem, 2)));
418 :
419 3 : if (auxX_Mem > ONEOVERS2C2CUTOFF[0]) {
420 0 : auxX_Mem = Lx * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
421 : } else {
422 3 : auxX_Mem = Lx * auxX_Mem / (2 * M_PI);
423 : }
424 :
425 3 : if (auxY_Mem > ONEOVERS2C2CUTOFF[0]) {
426 0 : auxY_Mem = Ly * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
427 : } else {
428 3 : auxY_Mem = Ly * auxY_Mem / (2 * M_PI);
429 : }
430 :
431 : // Number of lipid tails within the slice s of the membranes cylinder.
432 3 : std::vector<double> Nsp_Mem(NSMEM[0]), psi_Mem(NSMEM[0]), d_psi_Mem(NSMEM[0]);
433 :
434 : // Eq. 3 Hub & Awasthi JCTC 2017.
435 3 : double b_Mem = (ZETAMEM[0] / (1.0 - ZETAMEM[0])), c_Mem = ((1.0 - ZETAMEM[0]) * exp(b_Mem));
436 :
437 : // Eq. 19 Hub & Awasthi JCTC 2017.
438 3 : std::vector<double> fradial_Mem_d_faxial_Mem_dz(TAILS.size() * NSMEM[0]);
439 :
440 : // Eq. 20 Hub & Awasthi JCTC 2017.
441 3 : std::vector<double> Axs_Mem(NSMEM[0]), Ays_Mem(NSMEM[0]);
442 :
443 : // Eq. 1 Hub & Awasthi JCTC 2017. This is the CV that describes de Pore Nucleation.
444 3 : double Xi_Mem = 0.0;
445 :
446 : #ifdef _OPENMP
447 : #if _OPENMP >= 201307
448 3 : #pragma omp parallel for private(TailPosition,d_Xsc_Mem_dx,d_Xcc_Mem_dx,d_Ysc_Mem_dy,d_Ycc_Mem_dy,d_Xsc_Mem_dz,d_Xcc_Mem_dz,d_Ysc_Mem_dz,d_Ycc_Mem_dz,d_sx_Mem_dx,d_sy_Mem_dy,d_cx_Mem_dx,d_cy_Mem_dy,d_sx_Mem_dz,d_sy_Mem_dz,d_cx_Mem_dz,d_cy_Mem_dz,d_ws_Mem_dz,ri_Mem,x,fradial_Mem) reduction(vec_double_plus: Nsp_Mem, Axs_Mem, Ays_Mem)
449 : #endif
450 : #endif
451 : for (unsigned i = 0; i < TAILS.size(); i++) {
452 : if (analyzeThisParticle_Mem[i]) {
453 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
454 : d_Xsc_Mem_dx = 0.0;
455 : d_Xcc_Mem_dx = 0.0;
456 : d_Ysc_Mem_dy = 0.0;
457 : d_Ycc_Mem_dy = 0.0;
458 : d_Xsc_Mem_dz = 0.0;
459 : d_Xcc_Mem_dz = 0.0;
460 : d_Ysc_Mem_dz = 0.0;
461 : d_Ycc_Mem_dz = 0.0;
462 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
463 : if (Fs_Mem[s] != 0.0) {
464 : d_sx_Mem_dx = faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * cos(2.0 * M_PI * TailPosition[0]) / (Lx * Fs_Mem[s]);
465 : d_sy_Mem_dy = faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * cos(2.0 * M_PI * TailPosition[1]) / (Ly * Fs_Mem[s]);
466 : d_cx_Mem_dx = -faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * sin(2.0 * M_PI * TailPosition[0]) / (Lx * Fs_Mem[s]);
467 : d_cy_Mem_dy = -faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * sin(2.0 * M_PI * TailPosition[1]) / (Ly * Fs_Mem[s]);
468 : d_Xsc_Mem_dx += ws_Mem[s] * d_sx_Mem_dx / W_Mem;
469 : d_Xcc_Mem_dx += ws_Mem[s] * d_cx_Mem_dx / W_Mem;
470 : d_Ysc_Mem_dy += ws_Mem[s] * d_sy_Mem_dy / W_Mem;
471 : d_Ycc_Mem_dy += ws_Mem[s] * d_cy_Mem_dy / W_Mem;
472 :
473 : d_sx_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (sin(2.0 * M_PI * TailPosition[0]) - sx_Mem[s]) / Fs_Mem[s];
474 : d_sy_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (sin(2.0 * M_PI * TailPosition[1]) - sy_Mem[s]) / Fs_Mem[s];
475 : d_cx_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (cos(2.0 * M_PI * TailPosition[0]) - cx_Mem[s]) / Fs_Mem[s];
476 : d_cy_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (cos(2.0 * M_PI * TailPosition[1]) - cy_Mem[s]) / Fs_Mem[s];
477 : d_ws_Mem_dz = (1 - pow(ws_Mem[s], 2)) * d_faxial_Mem_dz[i + TAILS.size() * s];
478 : d_Xsc_Mem_dz += (ws_Mem[s] * d_sx_Mem_dz + d_ws_Mem_dz * (sx_Mem[s] - Xsc_Mem)) / W_Mem;
479 : d_Xcc_Mem_dz += (ws_Mem[s] * d_cx_Mem_dz + d_ws_Mem_dz * (cx_Mem[s] - Xcc_Mem)) / W_Mem;
480 : d_Ysc_Mem_dz += (ws_Mem[s] * d_sy_Mem_dz + d_ws_Mem_dz * (sy_Mem[s] - Ysc_Mem)) / W_Mem;
481 : d_Ycc_Mem_dz += (ws_Mem[s] * d_cy_Mem_dz + d_ws_Mem_dz * (cy_Mem[s] - Ycc_Mem)) / W_Mem;
482 : }
483 : }
484 : d_Xcyl_Mem_dx[i] = auxX_Mem * (-Xsc_Mem * d_Xcc_Mem_dx + Xcc_Mem * d_Xsc_Mem_dx);
485 : d_Xcyl_Mem_dz[i] = auxX_Mem * (-Xsc_Mem * d_Xcc_Mem_dz + Xcc_Mem * d_Xsc_Mem_dz);
486 : d_Ycyl_Mem_dy[i] = auxY_Mem * (-Ysc_Mem * d_Ycc_Mem_dy + Ycc_Mem * d_Ysc_Mem_dy);
487 : d_Ycyl_Mem_dz[i] = auxY_Mem * (-Ysc_Mem * d_Ycc_Mem_dz + Ycc_Mem * d_Ysc_Mem_dz);
488 :
489 : CylDistances_Mem[i] = pbcDistance(xyzCyl_Mem, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
490 : ri_Mem = sqrt(pow(CylDistances_Mem[i][0], 2) + pow(CylDistances_Mem[i][1], 2));
491 : x = ri_Mem / RCYLMEM[0];
492 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0]))) {
493 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0]))) {
494 : fradial_Mem = 1.0;
495 : } else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0]))) {
496 : fradial_Mem = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
497 : d_fradial_Mem_dx[i] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][0] / (RCYLMEM[0] * ri_Mem);
498 : d_fradial_Mem_dy[i] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][1] / (RCYLMEM[0] * ri_Mem);
499 : } else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0]))) {
500 : fradial_Mem = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
501 : d_fradial_Mem_dx[i] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][0] / (RCYLMEM[0] * ri_Mem);
502 : d_fradial_Mem_dy[i] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][1] / (RCYLMEM[0] * ri_Mem);
503 : }
504 :
505 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
506 : Nsp_Mem[s] += fradial_Mem * faxial_Mem[i + TAILS.size() * s];
507 : Axs_Mem[s] += faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dx[i];
508 : Ays_Mem[s] += faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dy[i];
509 : fradial_Mem_d_faxial_Mem_dz[i + TAILS.size() * s] = fradial_Mem * d_faxial_Mem_dz[i + TAILS.size() * s];
510 : }
511 : }
512 : }
513 : }
514 :
515 213 : for (unsigned s = 0; s < NSMEM[0]; s++) {
516 210 : if (Nsp_Mem[s] <= 1.0) {
517 166 : psi_Mem[s] = ZETAMEM[0] * Nsp_Mem[s];
518 166 : d_psi_Mem[s] = ZETAMEM[0];
519 166 : Xi_Mem += psi_Mem[s];
520 : } else {
521 44 : psi_Mem[s] = 1.0 - c_Mem * exp(-b_Mem * Nsp_Mem[s]);
522 44 : d_psi_Mem[s] = b_Mem * c_Mem * exp(-b_Mem * Nsp_Mem[s]);
523 44 : Xi_Mem += psi_Mem[s];
524 : }
525 : }
526 :
527 3 : Xi_Mem = Xi_Mem / NSMEM[0];
528 :
529 : // Eq. 18 Hub & Awasthi JCTC 2017.
530 3 : std::vector<double> faxial_Mem_d_fradial_Mem_dx(TAILS.size() * NSMEM[0]), faxial_Mem_d_fradial_Mem_dy(TAILS.size() * NSMEM[0]), faxial_Mem_d_fradial_Mem_dz(TAILS.size() * NSMEM[0]);
531 :
532 : // Eq. 13 Hub & Awasthi JCTC 2017.
533 3 : std::vector<Vector> derivatives_Mem(TAILS.size());
534 :
535 : #ifdef _OPENMP
536 : #if _OPENMP >= 201307
537 3 : #pragma omp parallel for private(aux)
538 : #endif
539 : #endif
540 : for (unsigned i = 0; i < TAILS.size(); i++) {
541 : if (analyzeThisParticle_Mem[i]) {
542 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
543 : if (faxial_Mem[i + TAILS.size() * s]) {
544 : faxial_Mem_d_fradial_Mem_dx[i + TAILS.size() * s] = faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dx[i] - d_Xcyl_Mem_dx[i] * Axs_Mem[s];
545 : faxial_Mem_d_fradial_Mem_dy[i + TAILS.size() * s] = faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dy[i] - d_Ycyl_Mem_dy[i] * Ays_Mem[s];
546 : faxial_Mem_d_fradial_Mem_dz[i + TAILS.size() * s] = -d_Xcyl_Mem_dz[i] * Axs_Mem[s] - d_Ycyl_Mem_dz[i] * Ays_Mem[s];
547 : }
548 : }
549 :
550 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
551 : aux = d_psi_Mem[s] / NSMEM[0];
552 : derivatives_Mem[i][0] += aux * faxial_Mem_d_fradial_Mem_dx[i + TAILS.size() * s];
553 : derivatives_Mem[i][1] += aux * faxial_Mem_d_fradial_Mem_dy[i + TAILS.size() * s];
554 : derivatives_Mem[i][2] += aux * (faxial_Mem_d_fradial_Mem_dz[i + TAILS.size() * s] + fradial_Mem_d_faxial_Mem_dz[i + TAILS.size() * s]);
555 : }
556 : }
557 : }
558 :
559 : // Derivatives and virial for the Xi_Mem.
560 3 : Tensor virial;
561 12291 : for (unsigned i = 0; i < TAILS.size(); i++) {
562 12288 : setAtomsDerivatives((i + membraneBeads), derivatives_Mem[i]);
563 12288 : virial -= Tensor(CylDistances_Mem[i], derivatives_Mem[i]);
564 : }
565 :
566 3 : setValue(Xi_Mem);
567 3 : setBoxDerivatives(virial);
568 3 : }
569 : }
570 : }
|
