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