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 FUSIONPOREEXPANSIONP
35 : /*
36 : A CV for inducing the expansion of a fusion pore from a nucleated fusion pore.
37 :
38 : Calculate the collective variable designed by Hub \cite Hub2021 and implemented into PLUMED by Masone and collaborators.
39 : This CV is capable of inducing the expansion of the fusion pore from a nucleated fusion pore.
40 :
41 : \f[
42 : \xi_e = \frac{R(r) - R_0}{R_0}
43 : \f]
44 :
45 : Where \f$\xi_e\f$ is the CV, \f$R_0\f$ is a normalization constant that makes zero the initial value of \f$\xi_e\f$, and
46 : \f$R(r)\f$ is the approximate radius of the fusion pore, which is defined by the number of waters and phosphateoxygens
47 : beads within a horizontal layer in the center of both membranes.
48 :
49 : \par Examples
50 :
51 : This example induces the expansion of a nucleated fusion pore (\f$\xi_e = 0.75\f$) from a just nucleated fusion pore (\f$\xi_e = 0.00\f$).
52 :
53 : \plumedfile
54 : lMem: GROUP ATOMS=1-10752,21505-22728,23953-24420 #All the lower membrane beads.
55 : uMem: GROUP ATOMS=10753-21504,22729-23952,24421-24888 #All the upper membrane beads.
56 : tails: GROUP ATOMS=8-23948:12,12-23952:12,23966-24884:18,23970-24888:18 #All the lipid tails beads (from the lower and upper membrane).
57 : waters: GROUP ATOMS=24889-56589 #All the water beads.
58 : po4: GROUP ATOMS=2-23942:12,23957-24875:18 #All the lipid phosphateoxygens beads.
59 :
60 : fusionPoreExpansion: FUSIONPOREEXPANSIONP UMEMBRANE=uMem LMEMBRANE=lMem TAILS=tails WATERS=waters PHOSPHATEOXYGENS=po4 NSMEM=85 D=7.0 R0=0.57
61 :
62 : MOVINGRESTRAINT ...
63 : ARG=fusionPoreExpansion
64 : STEP0=0 AT0=0.0 KAPPA0=10000.0
65 : STEP1=500000 AT1=0.75 KAPPA1=10000.0
66 : ...
67 :
68 : PRINT ARG=fusionPoreExpansion FILE=COLVAR STRIDE=1
69 :
70 : \endplumedfile
71 :
72 : */
73 : //+ENDPLUMEDOC
74 : class fusionPoreExpansionP : public Colvar
75 : {
76 : std::vector<AtomNumber> UMEM, LMEM, TAILS, WATERS, POXYGENS;
77 : std::vector<double> NSMEM, DSMEM, HMEM, VO, D, H, RMAX, R0, XCYL, YCYL;
78 :
79 : public:
80 : explicit fusionPoreExpansionP(const ActionOptions &);
81 : void calculate() override;
82 : static void registerKeywords(Keywords &keys);
83 : };
84 :
85 : PLUMED_REGISTER_ACTION(fusionPoreExpansionP, "FUSIONPOREEXPANSIONP")
86 :
87 3 : void fusionPoreExpansionP::registerKeywords(Keywords &keys)
88 : {
89 3 : Colvar::registerKeywords(keys);
90 6 : keys.add("atoms", "UMEMBRANE", "all the beads of the upper membrane.");
91 6 : keys.add("atoms", "LMEMBRANE", "all the beads of the lower membrane.");
92 6 : keys.add("atoms", "TAILS", "all the tail beads of the system.");
93 6 : keys.add("atoms", "WATERS", "all the water beads of the system.");
94 6 : keys.add("atoms", "PHOSPHATEOXYGENS", "all the lipid phosphateoxygens beads of the system.");
95 6 : keys.add("compulsory", "NSMEM", "the number of slices of the membrane fusion cylinder.");
96 6 : keys.add("optional", "DSMEM", "( default=0.1 ) thickness of the slices of the membrane fusion cylinder.");
97 6 : keys.add("optional", "HMEM", "( default=0.25 ) parameter of the step function θ(x,h) for the membrane fusion.");
98 6 : keys.add("optional", "VO", "( default=0.076879 ) beads' molecular volume.");
99 6 : keys.add("compulsory", "D", "horizontal layer thickness, it depends on the Z separation of the membranes.");
100 6 : keys.add("optional", "H", "( default=0.1 ) parameter of the step function θ(x,h) for the fusion pore expansion.");
101 6 : keys.add("optional", "RMAX", "( default=2.5 ) to avoid effects of membrane undulations in large membranes (more than 256 lipids).");
102 6 : keys.add("compulsory", "R0", "normalization constant that makes 0 the initial value of the CV.");
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", "X 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 : fusionPoreExpansionP::fusionPoreExpansionP(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 : parseAtomList("WATERS", WATERS);
123 1 : if (WATERS.size() == 0)
124 0 : error("WATERS has not any atom specified.");
125 :
126 2 : parseAtomList("PHOSPHATEOXYGENS", POXYGENS);
127 1 : if (POXYGENS.size() == 0)
128 0 : error("PHOSPHATEOXYGENS has not any atom specified.");
129 :
130 2 : parseVector("NSMEM", NSMEM);
131 1 : if (NSMEM.size() > 1)
132 0 : error("NSMEM cannot take more than one value.");
133 :
134 2 : parseVector("DSMEM", DSMEM);
135 1 : if (DSMEM.size() > 1)
136 0 : error("DSMEM cannot take more than one value.");
137 1 : if (DSMEM.size() == 0)
138 0 : DSMEM.push_back(0.1);
139 :
140 2 : parseVector("HMEM", HMEM);
141 1 : if (HMEM.size() > 1)
142 0 : error("HMEM cannot take more than one value.");
143 1 : if (HMEM.size() == 0)
144 0 : HMEM.push_back(0.25);
145 :
146 2 : parseVector("VO", VO);
147 1 : if (VO.size() > 1)
148 0 : error("VO cannot take more than one value.");
149 1 : if (VO.size() == 0)
150 0 : VO.push_back(0.076879);
151 :
152 2 : parseVector("D", D);
153 1 : if (D.size() > 1)
154 0 : error("D cannot take more than one value.");
155 :
156 2 : parseVector("H", H);
157 1 : if (H.size() > 1)
158 0 : error("H cannot take more than one value.");
159 1 : if (H.size() == 0)
160 0 : H.push_back(0.1);
161 :
162 2 : parseVector("RMAX", RMAX);
163 1 : if (RMAX.size() > 1)
164 0 : error("RMAX cannot take more than one value.");
165 1 : if (RMAX.size() == 0)
166 0 : RMAX.push_back(2.5);
167 :
168 2 : parseVector("R0", R0);
169 1 : if (R0.size() > 1)
170 0 : error("R0 cannot take more than one value.");
171 :
172 2 : parseVector("XCYL", XCYL);
173 1 : if (XCYL.size() > 1)
174 0 : error("XCYL cannot take more than one value.");
175 1 : if (XCYL.size() == 0)
176 1 : XCYL.push_back(-1.0);
177 :
178 2 : parseVector("YCYL", YCYL);
179 1 : if (YCYL.size() > 1)
180 0 : error("YCYL cannot take more than one value.");
181 1 : if (YCYL.size() == 0)
182 1 : YCYL.push_back(-1.0);
183 :
184 1 : checkRead();
185 :
186 : std::vector<AtomNumber> atoms;
187 12445 : for (unsigned i = 0; i < UMEM.size(); i++)
188 : {
189 12444 : atoms.push_back(UMEM[i]);
190 : }
191 12445 : for (unsigned i = 0; i < LMEM.size(); i++)
192 : {
193 12444 : atoms.push_back(LMEM[i]);
194 : }
195 4097 : for (unsigned i = 0; i < TAILS.size(); i++)
196 : {
197 4096 : atoms.push_back(TAILS[i]);
198 : }
199 31800 : for (unsigned i = 0; i < WATERS.size(); i++)
200 : {
201 31799 : atoms.push_back(WATERS[i]);
202 : }
203 2049 : for (unsigned i = 0; i < POXYGENS.size(); i++)
204 : {
205 2048 : atoms.push_back(POXYGENS[i]);
206 : }
207 :
208 1 : addValueWithDerivatives();
209 1 : setNotPeriodic();
210 1 : requestAtoms(atoms);
211 1 : }
212 4 : void fusionPoreExpansionP::calculate()
213 : {
214 : /*************************
215 : * *
216 : * System *
217 : * *
218 : **************************/
219 :
220 : // Box dimensions.
221 4 : double Lx = getBox()[0][0], Ly = getBox()[1][1], Lz = getBox()[2][2];
222 :
223 : // 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 .
224 : double ZuMem, ZuMemcos = 0.0, ZuMemsin = 0.0, uMemAngle, ZlMem, ZlMemcos = 0.0, ZlMemsin = 0.0, lMemAngle;
225 :
226 : #ifdef _OPENMP
227 : #if _OPENMP >= 201307
228 4 : #pragma omp parallel for private(uMemAngle, lMemAngle) reduction(+:ZuMemcos, ZuMemsin, ZlMemcos, ZlMemsin)
229 : #endif
230 : #endif
231 : for (unsigned i = 0; i < UMEM.size(); i++)
232 : {
233 : uMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
234 : lMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + UMEM.size())))[2];
235 : ZuMemcos += cos(uMemAngle);
236 : ZuMemsin += sin(uMemAngle);
237 : ZlMemcos += cos(lMemAngle);
238 : ZlMemsin += sin(lMemAngle);
239 : }
240 4 : ZuMemcos = ZuMemcos / UMEM.size();
241 4 : ZuMemsin = ZuMemsin / UMEM.size();
242 4 : ZuMem = Lz * (atan2(-ZuMemsin, -ZuMemcos) + M_PI) / (2.0 * M_PI);
243 4 : ZlMemcos = ZlMemcos / UMEM.size();
244 4 : ZlMemsin = ZlMemsin / UMEM.size();
245 4 : ZlMem = Lz * (atan2(-ZlMemsin, -ZlMemcos) + M_PI) / (2.0 * M_PI);
246 :
247 : // Z center of the boths membranes (upper and lower).
248 4 : double ZMems = (ZuMem + ZlMem) / 2.0;
249 :
250 : /**************************
251 : * *
252 : * Xcyl_Mem & Ycyl_Mem *
253 : * *
254 : ***************************/
255 :
256 : // Quantity of beads of the membranes.
257 4 : unsigned membraneBeads = UMEM.size() + LMEM.size();
258 :
259 : // Z distance from the lipid tail to the geometric center of both membranes.
260 : double ZTailDistance;
261 :
262 : // Z position of the first slice.
263 4 : double firstSliceZ_Mem = ZMems + (0.0 + 0.5 - NSMEM[0] / 2.0) * DSMEM[0];
264 :
265 : // Z distance between the first slice and the Z center of the membrane.
266 8 : double firstSliceZDist_Mem = pbcDistance(Vector(0.0, 0.0, firstSliceZ_Mem), Vector(0.0, 0.0, ZMems))[2];
267 :
268 : // Position in the cylinder.
269 : double PositionS_Mem;
270 :
271 : // Slices to analyze per particle.
272 : unsigned s1_Mem, s2_Mem;
273 :
274 : // Eq. 7 Hub & Awasthi JCTC 2017.
275 4 : std::vector<double> faxial_Mem(TAILS.size() * NSMEM[0]);
276 :
277 : // Eq. 10 Hub & Awasthi JCTC 2017.
278 4 : std::vector<double> Fs_Mem(NSMEM[0]);
279 :
280 : // Eq. 11 Hub & Awasthi JCTC 2017.
281 4 : std::vector<double> ws_Mem(NSMEM[0]);
282 :
283 : // Eq. 10 Hub & Awasthi JCTC 2017.
284 : double W_Mem = 0.0;
285 :
286 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
287 4 : std::vector<double> sx_Mem(NSMEM[0]), sy_Mem(NSMEM[0]), cx_Mem(NSMEM[0]), cy_Mem(NSMEM[0]);
288 :
289 : // Eq. 10 Hub & Awasthi JCTC 2017.
290 : double Xsc_Mem = 0.0, Xcc_Mem = 0.0, Ysc_Mem = 0.0, Ycc_Mem = 0.0;
291 :
292 : // Aux.
293 : double x, aux;
294 :
295 : // Scaled position of the lipid tail respect the origin of coordinates.
296 4 : Vector TailPosition;
297 :
298 : #ifdef _OPENMP
299 : #if _OPENMP >= 201307
300 : #pragma omp declare reduction(vec_double_plus : std::vector<double> : \
301 : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
302 : initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
303 : #endif
304 : #endif
305 :
306 : #ifdef _OPENMP
307 : #if _OPENMP >= 201307
308 4 : #pragma omp parallel for private(ZTailDistance, PositionS_Mem, TailPosition, x, aux, s1_Mem, s2_Mem) reduction(vec_double_plus:Fs_Mem, sx_Mem, sy_Mem, cx_Mem, cy_Mem)
309 : #endif
310 : #endif
311 : for (unsigned i = 0; i < TAILS.size(); i++)
312 : {
313 : ZTailDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + membraneBeads))[2];
314 : PositionS_Mem = (ZTailDistance + firstSliceZDist_Mem) / DSMEM[0];
315 : // If the following condition is met the particle is in the Z space of the cylinder.
316 : if ((PositionS_Mem >= (-0.5 - HMEM[0])) && (PositionS_Mem <= (NSMEM[0] + 0.5 - 1.0 + HMEM[0])))
317 : {
318 : //Defining the slices to analyze each particle.
319 : if (PositionS_Mem < 1)
320 : {
321 : s1_Mem = 0;
322 : s2_Mem = 2;
323 : }
324 : else if (PositionS_Mem <= (NSMEM[0] - 2.0))
325 : {
326 : s1_Mem = floor(PositionS_Mem) - 1;
327 : s2_Mem = floor(PositionS_Mem) + 1;
328 : }
329 : else
330 : {
331 : s1_Mem = NSMEM[0] - 3;
332 : s2_Mem = NSMEM[0] - 1;
333 : }
334 :
335 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
336 :
337 : for (unsigned s = s1_Mem; s <= s2_Mem; s++)
338 : {
339 : x = (ZTailDistance - (s + 0.5 - NSMEM[0] / 2.0) * DSMEM[0]) * 2.0 / DSMEM[0];
340 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0])))
341 : {
342 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0])))
343 : {
344 : faxial_Mem[i + TAILS.size() * s] = 1.0;
345 : Fs_Mem[s] += 1.0;
346 : sx_Mem[s] += sin(2.0 * M_PI * TailPosition[0]);
347 : sy_Mem[s] += sin(2.0 * M_PI * TailPosition[1]);
348 : cx_Mem[s] += cos(2.0 * M_PI * TailPosition[0]);
349 : cy_Mem[s] += cos(2.0 * M_PI * TailPosition[1]);
350 : }
351 : else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0])))
352 : {
353 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
354 : faxial_Mem[i + TAILS.size() * s] = aux;
355 : Fs_Mem[s] += aux;
356 : sx_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[0]);
357 : sy_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[1]);
358 : cx_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[0]);
359 : cy_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[1]);
360 : }
361 : else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0])))
362 : {
363 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
364 : faxial_Mem[i + TAILS.size() * s] = aux;
365 : Fs_Mem[s] += aux;
366 : sx_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[0]));
367 : sy_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[1]));
368 : cx_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[0]));
369 : cy_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[1]));
370 : }
371 : }
372 : }
373 : }
374 : }
375 :
376 344 : for (unsigned s = 0; s < NSMEM[0]; s++)
377 : {
378 340 : if (Fs_Mem[s] != 0.0)
379 : {
380 340 : ws_Mem[s] = tanh(Fs_Mem[s]);
381 340 : W_Mem += ws_Mem[s];
382 340 : sx_Mem[s] = sx_Mem[s] / Fs_Mem[s];
383 340 : sy_Mem[s] = sy_Mem[s] / Fs_Mem[s];
384 340 : cx_Mem[s] = cx_Mem[s] / Fs_Mem[s];
385 340 : cy_Mem[s] = cy_Mem[s] / Fs_Mem[s];
386 340 : Xsc_Mem += sx_Mem[s] * ws_Mem[s];
387 340 : Ysc_Mem += sy_Mem[s] * ws_Mem[s];
388 340 : Xcc_Mem += cx_Mem[s] * ws_Mem[s];
389 340 : Ycc_Mem += cy_Mem[s] * ws_Mem[s];
390 : }
391 : }
392 :
393 4 : Xsc_Mem = Xsc_Mem / W_Mem;
394 4 : Ysc_Mem = Ysc_Mem / W_Mem;
395 4 : Xcc_Mem = Xcc_Mem / W_Mem;
396 4 : Ycc_Mem = Ycc_Mem / W_Mem;
397 :
398 : // Eq. 12 Hub & Awasthi JCTC 2017.
399 : double Xcyl_Mem, Ycyl_Mem;
400 :
401 4 : if ((XCYL[0] > 0.0) && (YCYL[0] > 0.0))
402 : {
403 : Xcyl_Mem = XCYL[0];
404 : Ycyl_Mem = YCYL[0];
405 : }
406 : else
407 : {
408 4 : Xcyl_Mem = (atan2(-Xsc_Mem, -Xcc_Mem) + M_PI) * Lx / (2 * M_PI);
409 4 : Ycyl_Mem = (atan2(-Ysc_Mem, -Ycc_Mem) + M_PI) * Ly / (2 * M_PI);
410 : }
411 :
412 : /*************************
413 : * *
414 : * Xi_Exp *
415 : * *
416 : **************************/
417 :
418 : // Quantity of beads that could participate in the calculation of the Xi_Chain
419 4 : unsigned chainBeads = WATERS.size() + POXYGENS.size();
420 :
421 : // Quantity of beads that don't participate in the calculation of the Xi_Chain
422 4 : unsigned noChainBeads = (UMEM.size() * 2) + TAILS.size();
423 :
424 : // Center of the cylinder. X and Y are calculated (or defined), Z is the Z component of the geometric center of the membranes.
425 8 : Vector xyzCyl = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMems));
426 :
427 : // Estimation of RO with the Hub 2021 JCTC method. Only needed for the expansion.
428 4 : double RO = R0[0];
429 :
430 : // Number of polar atoms inside the horizontal layer. Eq. 3 Hub 2021 JCTC.
431 : double np = 0.0, fz, fr, fz_prime, fr_prime;
432 :
433 : // Derivative of np. Eq. 8 Hub 2021 JCTC.
434 4 : std::vector<double> d_np_dx(chainBeads), d_np_dy(chainBeads), d_np_dz(chainBeads);
435 :
436 : // Pore radius of the defect. Eq. 2 Hub 2021 JCTC.
437 : double poreR = 1.0;
438 :
439 : // Z center of the Membrane in the RMAX radius.
440 : double ZMemRMAX, ZMemRMAXcos = 0.0, ZMemRMAXsin = 0.0, countAux = 0.0, auxcos, auxsin;
441 :
442 : ZMemRMAX = ZMems;
443 :
444 : // The curvature of large membranes (1024 lipids) makes the Z-center of the membranes not to be representative
445 : // in some sectors, particularly in the region of the defect.
446 : //
447 : // To solve this, the center Z of the membranes in the defect sector is calculated and used to calculate
448 : // the number of polar atoms within the horizontal layer AND in the radious of the defect.
449 : //
450 : // ________ | | ________
451 : // ________ \_____| |______/ _______<-- Top membrane.
452 : // \______|P|_______/
453 : // |O|
454 : // | | <-- Z-center of the membranes in the region of the defect.
455 : // ______|R|_______ <-- Z-center of the membranes
456 : // / _____|E|______ \
457 : // / / | | \ \
458 : // ______/ / | | \ \______
459 : // _______/ \_______<-- Bottom membrane.
460 :
461 : // Center of mass for systems with PBC: https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions
462 4 : Vector MemCylDistances, distCylinder;
463 : double angle, ri;
464 :
465 : #ifdef _OPENMP
466 : #if _OPENMP >= 201307
467 4 : #pragma omp parallel for private(MemCylDistances, x, angle, auxcos, auxsin) reduction(+:ZMemRMAXcos, ZMemRMAXsin, countAux)
468 : #endif
469 : #endif
470 : for (unsigned i = 0; i < membraneBeads; i++)
471 : {
472 : MemCylDistances = pbcDistance(xyzCyl, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)));
473 : x = sqrt(pow(MemCylDistances[0], 2) + pow(MemCylDistances[1], 2)) / RMAX[0];
474 : if (!((x <= -1.0 - H[0]) || (x >= 1.0 + H[0])))
475 : {
476 : angle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
477 : auxcos = cos(angle);
478 : auxsin = sin(angle);
479 : if (((-1.0 + H[0]) <= x) && (x <= (1.0 - H[0])))
480 : {
481 : ZMemRMAXcos += 1.0 * auxcos;
482 : ZMemRMAXsin += 1.0 * auxsin;
483 : countAux += 1.0;
484 : }
485 : else if (((1.0 - H[0]) < x) && (x < (1.0 + H[0])))
486 : {
487 : ZMemRMAXcos += (0.5 - 0.75 * (x - 1.0) / H[0] + 0.25 * pow((x - 1.0), 3) / pow(H[0], 3)) * auxcos;
488 : ZMemRMAXsin += (0.5 - 0.75 * (x - 1.0) / H[0] + 0.25 * pow((x - 1.0), 3) / pow(H[0], 3)) * auxsin;
489 : countAux += (0.5 - 0.75 * (x - 1.0) / H[0] + 0.25 * pow((x - 1.0), 3) / pow(H[0], 3));
490 : }
491 : else if (((-1.0 - H[0]) < x) && (x < (-1.0 + H[0])))
492 : {
493 : ZMemRMAXcos += (0.5 + 0.75 * (x + 1.0) / H[0] - 0.25 * pow((x + 1.0), 3) / pow(H[0], 3)) * auxcos;
494 : ZMemRMAXsin += (0.5 + 0.75 * (x + 1.0) / H[0] - 0.25 * pow((x + 1.0), 3) / pow(H[0], 3)) * auxsin;
495 : countAux += (0.5 + 0.75 * (x + 1.0) / H[0] - 0.25 * pow((x + 1.0), 3) / pow(H[0], 3));
496 : }
497 : }
498 : }
499 :
500 4 : ZMemRMAXcos = ZMemRMAXcos / countAux;
501 4 : ZMemRMAXsin = ZMemRMAXsin / countAux;
502 4 : ZMemRMAX = Lz * (atan2(-ZMemRMAXsin, -ZMemRMAXcos) + M_PI) / (2.0 * M_PI);
503 :
504 4 : xyzCyl = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMemRMAX));
505 :
506 : #ifdef _OPENMP
507 : #if _OPENMP >= 201307
508 4 : #pragma omp parallel for private(distCylinder, fz, fz_prime, fr, fr_prime, ri, x) reduction(+:np)
509 : #endif
510 : #endif
511 : for (unsigned i = 0; i < chainBeads; i++)
512 : {
513 : distCylinder = pbcDistance(xyzCyl, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
514 : fz = 0.0;
515 : fz_prime = 0.0;
516 : fr = 0.0;
517 : fr_prime = 0.0;
518 :
519 : ri = sqrt(pow(distCylinder[0], 2) + pow(distCylinder[1], 2));
520 : x = ri / RMAX[0];
521 : if (!((x <= -1.0 - H[0]) || (x >= 1.0 + H[0])))
522 : {
523 : if (((-1.0 + H[0]) <= x) && (x <= (1.0 - H[0])))
524 : {
525 : fr = 1.0;
526 : }
527 : else if (((1.0 - H[0]) < x) && (x < (1.0 + H[0])))
528 : {
529 : fr = 0.5 - 0.75 * (x - 1.0) / H[0] + 0.25 * pow((x - 1.0), 3) / pow(H[0], 3);
530 : fr_prime = (-0.75 / H[0] + 0.75 * pow((x - 1.0), 2) / pow(H[0], 3)) / (RMAX[0] * ri);
531 : }
532 : else if (((-1.0 - H[0]) < x) && (x < (-1.0 + H[0])))
533 : {
534 : fr = 0.5 + 0.75 * (x + 1.0) / H[0] - 0.25 * pow((x + 1.0), 3) / pow(H[0], 3);
535 : fr_prime = (0.75 / H[0] - 0.75 * pow((x + 1), 2) / pow(H[0], 3)) / (RMAX[0] * ri);
536 : }
537 :
538 : x = distCylinder[2] * 2.0 / D[0];
539 : if (!((x <= -1.0 - H[0]) || (x >= 1.0 + H[0])))
540 : {
541 : if (((-1.0 + H[0]) <= x) && (x <= (1.0 - H[0])))
542 : {
543 : fz = 1.0;
544 : }
545 : else if (((1.0 - H[0]) < x) && (x < (1.0 + H[0])))
546 : {
547 : fz = 0.5 - 0.75 * (x - 1.0) / H[0] + 0.25 * pow((x - 1.0), 3) / pow(H[0], 3);
548 : fz_prime = (-0.75 / H[0] + 0.75 * pow((x - 1.0), 2) / pow(H[0], 3)) * 2.0 / D[0];
549 : }
550 : else if (((-1.0 - H[0]) < x) && (x < (-1.0 + H[0])))
551 : {
552 : fz = 0.5 + 0.75 * (x + 1.0) / H[0] - 0.25 * pow((x + 1.0), 3) / pow(H[0], 3);
553 : fz_prime = (0.75 / H[0] - 0.75 * pow((x + 1), 2) / pow(H[0], 3)) * 2.0 / D[0];
554 : }
555 :
556 : np += fz * fr;
557 : d_np_dx[i] = fz * fr_prime * distCylinder[0];
558 : d_np_dy[i] = fz * fr_prime * distCylinder[1];
559 : d_np_dz[i] = fz_prime * fr;
560 : }
561 : }
562 : }
563 4 : poreR = sqrt(np * VO[0] / (M_PI * D[0]));
564 :
565 : // This is the CV that describes the Pore Expansion.
566 4 : double Xi_Exp = (poreR - RO) / RO;
567 :
568 : // Derivatives vector.
569 4 : std::vector<Vector> derivatives(chainBeads);
570 :
571 : // Aux for the derivatives calculations. Eq. 7 Hub 2021 JCTC.
572 : double fact2 = 0.0;
573 :
574 4 : if (poreR != 0.0)
575 : {
576 4 : fact2 = VO[0] / (2.0 * M_PI * RO * D[0] * poreR);
577 : }
578 :
579 : // Distances from the oxygens to center of the cylinder.
580 4 : std::vector<Vector> CylDistances(chainBeads);
581 :
582 : #ifdef _OPENMP
583 : #if _OPENMP >= 201307
584 4 : #pragma omp parallel for
585 : #endif
586 : #endif
587 : for (unsigned i = 0; i < chainBeads; i++)
588 : {
589 : derivatives[i][0] = fact2 * d_np_dx[i];
590 : derivatives[i][1] = fact2 * d_np_dy[i];
591 : derivatives[i][2] = fact2 * d_np_dz[i];
592 : CylDistances[i] = pbcDistance(xyzCyl, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
593 : }
594 :
595 4 : Tensor virial;
596 135392 : for (unsigned i = 0; i < chainBeads; i++)
597 : {
598 135388 : setAtomsDerivatives((i + noChainBeads), derivatives[i]);
599 135388 : virial -= Tensor(CylDistances[i], derivatives[i]);
600 : }
601 :
602 4 : setValue(Xi_Exp);
603 4 : setBoxDerivatives(virial);
604 4 : }
605 : }
606 : }
|
