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