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