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