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 FUSIONPORENUCLEATIONP
37 : /*
38 : A CV for inducing the nucleation of the fusion pore from a hemifusion stalk.
39 :
40 : Calculate the collective variable designed by Hub and collaborators (see paper below) and
41 : implemented into PLUMED by Masone and collaborators.
42 : This CV is capable of inducing the nucleation of the fusion pore from a hemifusion stalk.
43 :
44 : $$
45 : \xi_n = \frac{1}{N_{sn}} \sum_{s=0}^{N_{sn}-1} \delta_{sn} (N_{sn}^{(p)})
46 : $$
47 :
48 : Where $\xi_n$ is the CV, $N_{sn}$ is the number of slices of the cylinder that make up the CV,
49 : $\delta_{sn}$ is a continuos function in the interval [0 1] ($\delta_{sf} = 0$ for no beads in the slice s, and
50 : $\delta_{sf} = 1$ for 1 or more beads in the slice s) and $N_{sf}^{(p)}$ accounts for the number of water and
51 : phosphateoxygens beads within the slice s.
52 :
53 : ## Examples
54 :
55 : This example induces the nucleation of the fusion pore ($\xi_n = 1.0$) from a hemifusion stalk ($\xi_n = 0.2$).
56 :
57 : ```plumed
58 :
59 : lMem: GROUP ATOMS=1-10752,21505-22728,23953-24420 #All the lower membrane beads.
60 : uMem: GROUP ATOMS=10753-21504,22729-23952,24421-24888 #All the upper membrane beads.
61 : 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).
62 : waters: GROUP ATOMS=24889-56490 #All the water beads.
63 : po4: GROUP ATOMS=2-23942:12,23957-24875:18 #All the lipid phosphateoxygens beads.
64 :
65 : fusionPoreNucleation: FUSIONPORENUCLEATIONP UMEMBRANE=uMem LMEMBRANE=lMem TAILS=tails WATERS=waters PHOSPHATEOXYGENS=po4 NSMEM=85 NS=45
66 :
67 : MOVINGRESTRAINT ...
68 : ARG=fusionPoreNucleation
69 : STEP0=0 AT0=0.2 KAPPA0=10000.0
70 : STEP1=500000 AT1=1.0 KAPPA1=10000.0
71 : ...
72 :
73 : PRINT ARG=fusionPoreNucleation FILE=COLVAR STRIDE=1
74 :
75 : ```
76 :
77 : */
78 : //+ENDPLUMEDOC
79 : class fusionPoreNucleationP : public Colvar {
80 : std::vector<AtomNumber> UMEM, LMEM, TAILS, WATERS, POXYGENS;
81 : std::vector<double> NSMEM, DSMEM, HMEM, NS, DS, HCH, RCYL, ZETA, ONEOVERS2C2CUTOFF, XCYL, YCYL;
82 :
83 : public:
84 : explicit fusionPoreNucleationP(const ActionOptions &);
85 : void calculate() override;
86 : static void registerKeywords(Keywords &keys);
87 : };
88 :
89 : PLUMED_REGISTER_ACTION(fusionPoreNucleationP, "FUSIONPORENUCLEATIONP")
90 :
91 3 : void fusionPoreNucleationP::registerKeywords(Keywords &keys) {
92 3 : Colvar::registerKeywords(keys);
93 3 : keys.add("atoms", "UMEMBRANE", "all the beads of the upper membrane.");
94 3 : keys.add("atoms", "LMEMBRANE", "all the beads of the lower membrane.");
95 3 : keys.add("atoms", "TAILS", "all the tail beads of the system.");
96 3 : keys.add("atoms", "WATERS", "all the water beads of the system.");
97 3 : keys.add("atoms", "PHOSPHATEOXYGENS", "all the lipid phosphateoxygens beads of the system.");
98 3 : keys.add("compulsory", "NSMEM", "the number of slices of the membrane fusion cylinder.");
99 3 : keys.add("optional", "DSMEM", "( default=0.1 ) thickness of the slices of the membrane fusion cylinder.");
100 3 : keys.add("optional", "HMEM", "( default=0.25 ) parameter of the step function θ(x,h) for the membrane fusion.");
101 3 : keys.add("compulsory", "NS", "the number of slices of the membrane-spanning cylinder in such a way that when the bilayers are flat and parallel the CV is equal to 0.2.");
102 3 : keys.add("optional", "DS", "( default=0.25 ) thickness of the slices of the membrane-spanning cylinder.");
103 3 : keys.add("optional", "HCH", "( default=0.25 ) parameter of the step function θ(x,h) for the CV.");
104 3 : keys.add("optional", "RCYL", "( default=0.8 ) the radius of the membrane-spanning cylinder.");
105 3 : keys.add("optional", "ZETA", "( default=0.75 ) parameter of the switch function ψ(x,ζ).");
106 3 : keys.add("optional", "ONEOVERS2C2CUTOFF", "( default=500 ) cut off large values for the derivative of the atan2 function to avoid violate energy.");
107 3 : keys.add("optional", "XCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
108 3 : keys.add("optional", "YCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
109 6 : keys.setValueDescription("scalar","the value of the CV");
110 3 : keys.addDOI("10.1021/acs.jctc.7b00106");
111 3 : }
112 :
113 1 : fusionPoreNucleationP::fusionPoreNucleationP(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao) {
114 2 : parseAtomList("UMEMBRANE", UMEM);
115 1 : if (UMEM.size() == 0) {
116 0 : error("UMEMBRANE has not any atom specified.");
117 : }
118 :
119 2 : parseAtomList("LMEMBRANE", LMEM);
120 1 : if (LMEM.size() == 0) {
121 0 : error("LMEMBRANE has not any atom specified.");
122 : }
123 :
124 2 : parseAtomList("TAILS", TAILS);
125 1 : if (TAILS.size() == 0) {
126 0 : error("TAILS has not any atom specified.");
127 : }
128 :
129 2 : parseAtomList("WATERS", WATERS);
130 1 : if (WATERS.size() == 0) {
131 0 : error("WATERS has not any atom specified.");
132 : }
133 :
134 2 : parseAtomList("PHOSPHATEOXYGENS", POXYGENS);
135 1 : if (POXYGENS.size() == 0) {
136 0 : error("PHOSPHATEOXYGENS has not any atom specified.");
137 : }
138 :
139 2 : parseVector("NSMEM", NSMEM);
140 1 : if (NSMEM.size() > 1) {
141 0 : error("NSMEM cannot take more than one value.");
142 : }
143 :
144 2 : parseVector("DSMEM", DSMEM);
145 1 : if (DSMEM.size() > 1) {
146 0 : error("DSMEM cannot take more than one value.");
147 : }
148 1 : if (DSMEM.size() == 0) {
149 0 : DSMEM.push_back(0.1);
150 : }
151 :
152 2 : parseVector("HMEM", HMEM);
153 1 : if (HMEM.size() > 1) {
154 0 : error("HMEM cannot take more than one value.");
155 : }
156 1 : if (HMEM.size() == 0) {
157 0 : HMEM.push_back(0.25);
158 : }
159 :
160 2 : parseVector("NS", NS);
161 1 : if (NS.size() > 1) {
162 0 : error("NS cannot take more than one value.");
163 : }
164 :
165 2 : parseVector("DS", DS);
166 1 : if (DS.size() > 1) {
167 0 : error("DS cannot take more than one value.");
168 : }
169 1 : if (DS.size() == 0) {
170 0 : DS.push_back(0.25);
171 : }
172 :
173 2 : parseVector("HCH", HCH);
174 1 : if (HCH.size() > 1) {
175 0 : error("H cannot take more than one value.");
176 : }
177 1 : if (HCH.size() == 0) {
178 0 : HCH.push_back(0.25);
179 : }
180 :
181 2 : parseVector("RCYL", RCYL);
182 1 : if (RCYL.size() > 1) {
183 0 : error("RCYL cannot take more than one value.");
184 : }
185 1 : if (RCYL.size() == 0) {
186 0 : RCYL.push_back(0.8);
187 : }
188 :
189 2 : parseVector("ZETA", ZETA);
190 1 : if (ZETA.size() > 1) {
191 0 : error("ZETA cannot take more than one value.");
192 : }
193 1 : if (ZETA.size() == 0) {
194 0 : ZETA.push_back(0.75);
195 : }
196 :
197 2 : parseVector("ONEOVERS2C2CUTOFF", ONEOVERS2C2CUTOFF);
198 1 : if (ONEOVERS2C2CUTOFF.size() > 1) {
199 0 : error("ONEOVERS2C2CUTOFF cannot take more than one value.");
200 : }
201 1 : if (ONEOVERS2C2CUTOFF.size() == 0) {
202 1 : ONEOVERS2C2CUTOFF.push_back(500);
203 : }
204 :
205 2 : parseVector("XCYL", XCYL);
206 1 : if (XCYL.size() > 1) {
207 0 : error("XCYL cannot take more than one value.");
208 : }
209 1 : if (XCYL.size() == 0) {
210 1 : XCYL.push_back(-1.0);
211 : }
212 :
213 2 : parseVector("YCYL", YCYL);
214 1 : if (YCYL.size() > 1) {
215 0 : error("YCYL cannot take more than one value.");
216 : }
217 1 : if (YCYL.size() == 0) {
218 1 : YCYL.push_back(-1.0);
219 : }
220 :
221 1 : checkRead();
222 :
223 : std::vector<AtomNumber> atoms;
224 12445 : for (unsigned i = 0; i < UMEM.size(); i++) {
225 12444 : atoms.push_back(UMEM[i]);
226 : }
227 12445 : for (unsigned i = 0; i < LMEM.size(); i++) {
228 12444 : atoms.push_back(LMEM[i]);
229 : }
230 4097 : for (unsigned i = 0; i < TAILS.size(); i++) {
231 4096 : atoms.push_back(TAILS[i]);
232 : }
233 31603 : for (unsigned i = 0; i < WATERS.size(); i++) {
234 31602 : atoms.push_back(WATERS[i]);
235 : }
236 2049 : for (unsigned i = 0; i < POXYGENS.size(); i++) {
237 2048 : atoms.push_back(POXYGENS[i]);
238 : }
239 :
240 1 : addValueWithDerivatives();
241 1 : setNotPeriodic();
242 1 : requestAtoms(atoms);
243 1 : }
244 :
245 4 : void fusionPoreNucleationP::calculate() {
246 : /*************************
247 : * *
248 : * System *
249 : * *
250 : **************************/
251 :
252 : // Box dimensions.
253 4 : double Lx = getBox()[0][0], Ly = getBox()[1][1], Lz = getBox()[2][2];
254 :
255 : // 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 .
256 : double ZuMem, ZuMemcos = 0.0, ZuMemsin = 0.0, uMemAngle, ZlMem, ZlMemcos = 0.0, ZlMemsin = 0.0, lMemAngle;
257 :
258 : #ifdef _OPENMP
259 : #if _OPENMP >= 201307
260 4 : #pragma omp parallel for private(uMemAngle, lMemAngle) reduction(+:ZuMemcos, ZuMemsin, ZlMemcos, ZlMemsin)
261 : #endif
262 : #endif
263 : for (unsigned i = 0; i < UMEM.size(); i++) {
264 : uMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
265 : lMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + UMEM.size())))[2];
266 : ZuMemcos += cos(uMemAngle);
267 : ZuMemsin += sin(uMemAngle);
268 : ZlMemcos += cos(lMemAngle);
269 : ZlMemsin += sin(lMemAngle);
270 : }
271 4 : ZuMemcos = ZuMemcos / UMEM.size();
272 4 : ZuMemsin = ZuMemsin / UMEM.size();
273 4 : ZuMem = Lz * (atan2(-ZuMemsin, -ZuMemcos) + M_PI) / (2.0 * M_PI);
274 4 : ZlMemcos = ZlMemcos / UMEM.size();
275 4 : ZlMemsin = ZlMemsin / UMEM.size();
276 4 : ZlMem = Lz * (atan2(-ZlMemsin, -ZlMemcos) + M_PI) / (2.0 * M_PI);
277 :
278 : // Z center of the boths membranes (upper and lower).
279 4 : double ZMems = (ZuMem + ZlMem) / 2.0;
280 :
281 : /**************************
282 : * *
283 : * Xcyl_Mem & Ycyl_Mem *
284 : * *
285 : ***************************/
286 :
287 : // Quantity of beads of the membranes.
288 4 : unsigned membraneBeads = UMEM.size() + LMEM.size();
289 :
290 : // Z distance from the lipid tail to the geometric center of both membranes.
291 : double ZTailDistance;
292 :
293 : // Z position of the first slice.
294 4 : double firstSliceZ_Mem = ZMems + (0.0 + 0.5 - NSMEM[0] / 2.0) * DSMEM[0];
295 :
296 : // Z distance between the first slice and the Z center of the membrane.
297 4 : double firstSliceZDist_Mem = pbcDistance(Vector(0.0, 0.0, firstSliceZ_Mem), Vector(0.0, 0.0, ZMems))[2];
298 :
299 : // Position in the cylinder.
300 : double PositionS_Mem;
301 :
302 : // Slices to analyze per particle.
303 : unsigned s1_Mem, s2_Mem;
304 :
305 : // Eq. 7 Hub & Awasthi JCTC 2017.
306 4 : std::vector<double> faxial_Mem(TAILS.size() * NSMEM[0]);
307 :
308 : // Eq. 10 Hub & Awasthi JCTC 2017.
309 4 : std::vector<double> Fs_Mem(NSMEM[0]);
310 :
311 : // Eq. 11 Hub & Awasthi JCTC 2017.
312 4 : std::vector<double> ws_Mem(NSMEM[0]);
313 :
314 : // Eq. 10 Hub & Awasthi JCTC 2017.
315 : double W_Mem = 0.0;
316 :
317 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
318 4 : std::vector<double> sx_Mem(NSMEM[0]), sy_Mem(NSMEM[0]), cx_Mem(NSMEM[0]), cy_Mem(NSMEM[0]);
319 :
320 : // Eq. 10 Hub & Awasthi JCTC 2017.
321 : double Xsc_Mem = 0.0, Xcc_Mem = 0.0, Ysc_Mem = 0.0, Ycc_Mem = 0.0;
322 :
323 : // Aux.
324 : double x, aux;
325 :
326 : // Scaled position of the lipid tail respect the origin of coordinates.
327 : Vector TailPosition;
328 :
329 : #ifdef _OPENMP
330 : #if _OPENMP >= 201307
331 : #pragma omp declare reduction(vec_double_plus : std::vector<double> : \
332 : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
333 : initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
334 : #endif
335 : #endif
336 :
337 : #ifdef _OPENMP
338 : #if _OPENMP >= 201307
339 4 : #pragma omp parallel for private(ZTailDistance, PositionS_Mem, s1_Mem, s2_Mem, TailPosition, x, aux) reduction(vec_double_plus:Fs_Mem, sx_Mem, sy_Mem, cx_Mem, cy_Mem)
340 : #endif
341 : #endif
342 : for (unsigned i = 0; i < TAILS.size(); i++) {
343 : ZTailDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + membraneBeads))[2];
344 : PositionS_Mem = (ZTailDistance + firstSliceZDist_Mem) / DSMEM[0];
345 : // If the following condition is met the particle is in the Z space of the cylinder.
346 : if ((PositionS_Mem >= (-0.5 - HMEM[0])) && (PositionS_Mem <= (NSMEM[0] + 0.5 - 1.0 + HMEM[0]))) {
347 : //Defining the slices to analyze each particle.
348 : if (PositionS_Mem < 1) {
349 : s1_Mem = 0;
350 : s2_Mem = 2;
351 : } else if (PositionS_Mem <= (NSMEM[0] - 2.0)) {
352 : s1_Mem = floor(PositionS_Mem) - 1;
353 : s2_Mem = floor(PositionS_Mem) + 1;
354 : } else {
355 : s1_Mem = NSMEM[0] - 3;
356 : s2_Mem = NSMEM[0] - 1;
357 : }
358 :
359 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
360 :
361 : for (unsigned s = s1_Mem; s <= s2_Mem; s++) {
362 : x = (ZTailDistance - (s + 0.5 - NSMEM[0] / 2.0) * DSMEM[0]) * 2.0 / DSMEM[0];
363 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0]))) {
364 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0]))) {
365 : faxial_Mem[i + TAILS.size() * s] = 1.0;
366 : Fs_Mem[s] += 1.0;
367 : sx_Mem[s] += sin(2.0 * M_PI * TailPosition[0]);
368 : sy_Mem[s] += sin(2.0 * M_PI * TailPosition[1]);
369 : cx_Mem[s] += cos(2.0 * M_PI * TailPosition[0]);
370 : cy_Mem[s] += cos(2.0 * M_PI * TailPosition[1]);
371 : } else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0]))) {
372 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
373 : faxial_Mem[i + TAILS.size() * s] = aux;
374 : Fs_Mem[s] += aux;
375 : sx_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[0]);
376 : sy_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[1]);
377 : cx_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[0]);
378 : cy_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[1]);
379 : } else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0]))) {
380 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
381 : faxial_Mem[i + TAILS.size() * s] = aux;
382 : Fs_Mem[s] += aux;
383 : sx_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[0]));
384 : sy_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[1]));
385 : cx_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[0]));
386 : cy_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[1]));
387 : }
388 : }
389 : }
390 : }
391 : }
392 :
393 344 : for (unsigned s = 0; s < NSMEM[0]; s++) {
394 340 : if (Fs_Mem[s] != 0.0) {
395 339 : ws_Mem[s] = tanh(Fs_Mem[s]);
396 339 : W_Mem += ws_Mem[s];
397 339 : sx_Mem[s] = sx_Mem[s] / Fs_Mem[s];
398 339 : sy_Mem[s] = sy_Mem[s] / Fs_Mem[s];
399 339 : cx_Mem[s] = cx_Mem[s] / Fs_Mem[s];
400 339 : cy_Mem[s] = cy_Mem[s] / Fs_Mem[s];
401 339 : Xsc_Mem += sx_Mem[s] * ws_Mem[s];
402 339 : Ysc_Mem += sy_Mem[s] * ws_Mem[s];
403 339 : Xcc_Mem += cx_Mem[s] * ws_Mem[s];
404 339 : Ycc_Mem += cy_Mem[s] * ws_Mem[s];
405 : }
406 : }
407 :
408 4 : Xsc_Mem = Xsc_Mem / W_Mem;
409 4 : Ysc_Mem = Ysc_Mem / W_Mem;
410 4 : Xcc_Mem = Xcc_Mem / W_Mem;
411 4 : Ycc_Mem = Ycc_Mem / W_Mem;
412 :
413 : // Eq. 12 Hub & Awasthi JCTC 2017.
414 : double Xcyl_Mem, Ycyl_Mem;
415 :
416 4 : if ((XCYL[0] > 0.0) && (YCYL[0] > 0.0)) {
417 : Xcyl_Mem = XCYL[0];
418 : Ycyl_Mem = YCYL[0];
419 : } else {
420 4 : Xcyl_Mem = (atan2(-Xsc_Mem, -Xcc_Mem) + M_PI) * Lx / (2 * M_PI);
421 4 : Ycyl_Mem = (atan2(-Ysc_Mem, -Ycc_Mem) + M_PI) * Ly / (2 * M_PI);
422 : }
423 :
424 : /*************************
425 : * *
426 : * Xi_n *
427 : * *
428 : **************************/
429 :
430 : // Eq. 1 Hub & Awasthi JCTC 2017. This is the CV that describes de Pore Nucleation.
431 4 : double Xi_n = 0.0;
432 :
433 : // Quantity of beads that could participate in the calculation of the Xi_Chain
434 4 : unsigned chainBeads = WATERS.size() + POXYGENS.size();
435 :
436 : // Quantity of beads that don't participate in the calculation of the Xi_Chain
437 4 : unsigned noChainBeads = (UMEM.size() * 2) + TAILS.size();
438 :
439 : // Z Distances from the oxygens to the geometric center of the membranes.
440 : double ZMemDistance;
441 :
442 : // Scaled positions of the oxygens to respect of the origin of coordinates.
443 : Vector Position;
444 :
445 : // Distance from the water/phosphate group to the defect cylinder.
446 : Vector distCylinder;
447 :
448 : // Center of the cylinder. XY components are calculated (or defined), Z is the Z geometric center of the membranes of the system.
449 4 : Vector xyzCyl_Mem = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMems));
450 :
451 : // Average of the radius of the water and lipid cylinder.
452 4 : double RCYLAVERAGE = RCYL[0] * (1 + HCH[0]);
453 :
454 : // Conditions.
455 : bool condition1, condition2, condition3;
456 :
457 : // Z position of the first slice.
458 4 : double firstSliceZ = ZMems + (0.0 + 0.5 - NS[0] / 2.0) * DS[0];
459 :
460 : // Z distance between the first slice and the Z center of the membrane.
461 4 : double firstSliceZDist = pbcDistance(Vector(0.0, 0.0, firstSliceZ), Vector(0.0, 0.0, ZMems))[2];
462 :
463 : // Position in the cylinder.
464 : double PositionS;
465 :
466 : // Mark the particles to analyze.
467 4 : std::vector<double> analyzeThisParticle(chainBeads);
468 :
469 : // Slices to analyze per particle.
470 4 : std::vector<unsigned> s1(chainBeads), s2(chainBeads);
471 :
472 : // Eq. 7 Hub & Awasthi JCTC 2017.
473 4 : std::vector<double> faxial(chainBeads * NS[0]);
474 :
475 : // Eq. 16 Hub & Awasthi JCTC 2017.
476 4 : std::vector<double> d_faxial_dz(chainBeads * NS[0]);
477 :
478 : // Eq. 10 Hub & Awasthi JCTC 2017.
479 4 : std::vector<double> Fs(NS[0]);
480 :
481 : // Eq. 11 Hub & Awasthi JCTC 2017.
482 4 : std::vector<double> ws(NS[0]);
483 :
484 : // Eq. 10 Hub & Awasthi JCTC 2017.
485 : double W = 0.0;
486 :
487 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
488 4 : std::vector<double> sx(NS[0]), sy(NS[0]), cx(NS[0]), cy(NS[0]);
489 :
490 : // Eq. 10 Hub & Awasthi JCTC 2017.
491 : double Xsc = 0.0, Xcc = 0.0, Ysc = 0.0, Ycc = 0.0;
492 :
493 : #ifdef _OPENMP
494 : #if _OPENMP >= 201307
495 4 : #pragma omp parallel for private(distCylinder, aux, condition1, condition2, condition3, ZMemDistance, PositionS, Position, x) reduction(vec_double_plus:Fs, sx, sy, cx, cy)
496 : #endif
497 : #endif
498 : for (unsigned i = 0; i < chainBeads; i++) {
499 : distCylinder = pbcDistance(xyzCyl_Mem, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
500 : aux = sqrt(pow(distCylinder[0], 2) + pow(distCylinder[1], 2));
501 : condition1 = ((aux / RCYLAVERAGE) < 1.0);
502 : condition2 = ((pbcDistance(Vector(0.0, 0.0, ZuMem), getPosition(i + noChainBeads))[2] > 0) && (aux / RCYLAVERAGE) < 2.0);
503 : condition3 = ((pbcDistance(getPosition(i + noChainBeads), Vector(0.0, 0.0, ZlMem))[2] > 0) && (aux / RCYLAVERAGE) < 2.0);
504 : if (condition1 || condition2 || condition3) {
505 : ZMemDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + noChainBeads))[2];
506 : PositionS = (ZMemDistance + firstSliceZDist) / DS[0];
507 : // If the following condition is met the particle is in the Z space of the cylinder.
508 : if ((PositionS >= (-0.5 - HCH[0])) && (PositionS <= (NS[0] + 0.5 - 1.0 + HCH[0]))) {
509 : analyzeThisParticle[i] = 1.0;
510 :
511 : //Defining the slices to analyze each particle.
512 : if (PositionS < 1) {
513 : s1[i] = 0;
514 : s2[i] = 2;
515 : } else if (PositionS <= (NS[0] - 2.0)) {
516 : s1[i] = floor(PositionS) - 1;
517 : s2[i] = floor(PositionS) + 1;
518 : } else {
519 : s1[i] = NS[0] - 3;
520 : s2[i] = NS[0] - 1;
521 : }
522 :
523 : Position = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
524 :
525 : for (unsigned s = s1[i]; s <= s2[i]; s++) {
526 : x = (ZMemDistance - (s + 0.5 - NS[0] / 2.0) * DS[0]) * 2.0 / DS[0];
527 : if (!((x <= -1.0 - HCH[0]) || (x >= 1.0 + HCH[0]))) {
528 : if (((-1.0 + HCH[0]) <= x) && (x <= (1.0 - HCH[0]))) {
529 : faxial[i + chainBeads * s] = 1.0;
530 : Fs[s] += 1.0;
531 : sx[s] += sin(2.0 * M_PI * Position[0]);
532 : sy[s] += sin(2.0 * M_PI * Position[1]);
533 : cx[s] += cos(2.0 * M_PI * Position[0]);
534 : cy[s] += cos(2.0 * M_PI * Position[1]);
535 : } else if (((1.0 - HCH[0]) < x) && (x < (1.0 + HCH[0]))) {
536 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HCH[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HCH[0], 3)));
537 : faxial[i + chainBeads * s] = aux;
538 : d_faxial_dz[i + chainBeads * s] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * 2.0 / DS[0];
539 : Fs[s] += aux;
540 : sx[s] += aux * sin(2.0 * M_PI * Position[0]);
541 : sy[s] += aux * sin(2.0 * M_PI * Position[1]);
542 : cx[s] += aux * cos(2.0 * M_PI * Position[0]);
543 : cy[s] += aux * cos(2.0 * M_PI * Position[1]);
544 : } else if (((-1.0 - HCH[0]) < x) && (x < (-1.0 + HCH[0]))) {
545 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HCH[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HCH[0], 3)));
546 : faxial[i + chainBeads * s] = aux;
547 : d_faxial_dz[i + chainBeads * s] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * 2.0 / DS[0];
548 : Fs[s] += aux;
549 : sx[s] += (aux * sin(2.0 * M_PI * Position[0]));
550 : sy[s] += (aux * sin(2.0 * M_PI * Position[1]));
551 : cx[s] += (aux * cos(2.0 * M_PI * Position[0]));
552 : cy[s] += (aux * cos(2.0 * M_PI * Position[1]));
553 : }
554 : }
555 : }
556 : }
557 : }
558 : }
559 :
560 184 : for (unsigned s = 0; s < NS[0]; s++) {
561 180 : if (Fs[s] != 0.0) {
562 49 : ws[s] = tanh(Fs[s]);
563 49 : W += ws[s];
564 49 : sx[s] = sx[s] / Fs[s];
565 49 : sy[s] = sy[s] / Fs[s];
566 49 : cx[s] = cx[s] / Fs[s];
567 49 : cy[s] = cy[s] / Fs[s];
568 49 : Xsc += sx[s] * ws[s];
569 49 : Ysc += sy[s] * ws[s];
570 49 : Xcc += cx[s] * ws[s];
571 49 : Ycc += cy[s] * ws[s];
572 : }
573 : }
574 :
575 4 : Xsc = Xsc / W;
576 4 : Ysc = Ysc / W;
577 4 : Xcc = Xcc / W;
578 4 : Ycc = Ycc / W;
579 :
580 : // Eq. 12 Hub & Awasthi JCTC 2017.
581 : double Xcyl, Ycyl;
582 :
583 : Xcyl = Xcyl_Mem;
584 : Ycyl = Ycyl_Mem;
585 :
586 : // Eq. 25, 26 and 27 Hub & Awasthi JCTC 2017.
587 : double d_sx_dx, d_sx_dz, d_sy_dy, d_sy_dz, d_cx_dx, d_cx_dz, d_cy_dy, d_cy_dz;
588 :
589 : // Eq. 29 Hub & Awasthi JCTC 2017.
590 : double d_ws_dz;
591 :
592 : // Eq. 31, 32 and 33 Hub & Awasthi JCTC 2017
593 : double d_Xsc_dx, d_Xsc_dz, d_Xcc_dx, d_Xcc_dz, d_Ysc_dy, d_Ysc_dz, d_Ycc_dy, d_Ycc_dz;
594 :
595 : // Center of the cylinder. X and Y are calculated (or defined), Z is the Z component of the geometric center of the membranes.
596 4 : Vector xyzCyl = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl, Ycyl, ZMems));
597 :
598 : // Distances from the oxygens to center of the cylinder.
599 4 : std::vector<Vector> CylDistances(chainBeads);
600 :
601 : // Modulo of the XY distances from the oxygens to the center of the cylinder.
602 : double ri;
603 :
604 : // Eq. 8 Hub & Awasthi JCTC 2017.
605 : double fradial;
606 :
607 : // Eq. 15 Hub & Awasthi JCTC 2017.
608 4 : std::vector<double> d_fradial_dx(chainBeads), d_fradial_dy(chainBeads);
609 :
610 : // Eq. 35, 36, 37 and 38 Hub & Awasthi JCTC 2017.
611 4 : std::vector<double> d_Xcyl_dx(chainBeads), d_Xcyl_dz(chainBeads), d_Ycyl_dy(chainBeads), d_Ycyl_dz(chainBeads);
612 :
613 : // To avoid rare instabilities auxX and auxY are truncated at a configurable value (default 500).
614 4 : double auxX = (1 / (pow(Xsc, 2) + pow(Xcc, 2))), auxY = (1 / (pow(Ysc, 2) + pow(Ycc, 2)));
615 :
616 4 : if (auxX > ONEOVERS2C2CUTOFF[0]) {
617 0 : auxX = Lx * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
618 : } else {
619 4 : auxX = Lx * auxX / (2 * M_PI);
620 : }
621 :
622 4 : if (auxY > ONEOVERS2C2CUTOFF[0]) {
623 0 : auxY = Ly * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
624 : } else {
625 4 : auxY = Ly * auxY / (2 * M_PI);
626 : }
627 :
628 : //Number of oxygens within the slice s of the membrane-spanning cylinder.
629 4 : std::vector<double> Nsp(NS[0]), psi(NS[0]), d_psi(NS[0]);
630 :
631 : // Eq. 3 Hub & Awasthi JCTC 2017.
632 4 : double b = (ZETA[0] / (1.0 - ZETA[0])), c = ((1.0 - ZETA[0]) * exp(b));
633 :
634 : // Eq. 19 Hub & Awasthi JCTC 2017.
635 4 : std::vector<double> fradial_d_faxial_dz(chainBeads * NS[0]);
636 :
637 : // Eq. 20 Hub & Awasthi JCTC 2017.
638 4 : std::vector<double> Axs(NS[0]), Ays(NS[0]);
639 :
640 : #ifdef _OPENMP
641 : #if _OPENMP >= 201307
642 4 : #pragma omp parallel for private(d_Xsc_dx,d_Xcc_dx,d_Ysc_dy,d_Ycc_dy,d_Xsc_dz,d_Xcc_dz,d_Ysc_dz,d_Ycc_dz,d_sx_dx,d_sy_dy,d_cx_dx,d_cy_dy,d_sx_dz,d_sy_dz,d_cx_dz,d_cy_dz,d_ws_dz,ri,x,fradial) reduction(vec_double_plus: Nsp, Axs, Ays)
643 : #endif
644 : #endif
645 : for (unsigned i = 0; i < chainBeads; i++) {
646 : CylDistances[i] = pbcDistance(xyzCyl, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
647 : if (analyzeThisParticle[i]) {
648 : Position = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
649 : d_Xsc_dx = 0.0;
650 : d_Xcc_dx = 0.0;
651 : d_Ysc_dy = 0.0;
652 : d_Ycc_dy = 0.0;
653 : d_Xsc_dz = 0.0;
654 : d_Xcc_dz = 0.0;
655 : d_Ysc_dz = 0.0;
656 : d_Ycc_dz = 0.0;
657 : for (unsigned s = s1[i]; s <= s2[i]; s++) {
658 : if (Fs[s] != 0.0) {
659 : d_sx_dx = faxial[i + chainBeads * s] * 2.0 * M_PI * cos(2.0 * M_PI * Position[0]) / (Lx * Fs[s]);
660 : d_sy_dy = faxial[i + chainBeads * s] * 2.0 * M_PI * cos(2.0 * M_PI * Position[1]) / (Ly * Fs[s]);
661 : d_cx_dx = -faxial[i + chainBeads * s] * 2.0 * M_PI * sin(2.0 * M_PI * Position[0]) / (Lx * Fs[s]);
662 : d_cy_dy = -faxial[i + chainBeads * s] * 2.0 * M_PI * sin(2.0 * M_PI * Position[1]) / (Ly * Fs[s]);
663 : d_Xsc_dx += ws[s] * d_sx_dx / W;
664 : d_Xcc_dx += ws[s] * d_cx_dx / W;
665 : d_Ysc_dy += ws[s] * d_sy_dy / W;
666 : d_Ycc_dy += ws[s] * d_cy_dy / W;
667 :
668 : d_sx_dz = d_faxial_dz[i + chainBeads * s] * (sin(2.0 * M_PI * Position[0]) - sx[s]) / Fs[s];
669 : d_sy_dz = d_faxial_dz[i + chainBeads * s] * (sin(2.0 * M_PI * Position[1]) - sy[s]) / Fs[s];
670 : d_cx_dz = d_faxial_dz[i + chainBeads * s] * (cos(2.0 * M_PI * Position[0]) - cx[s]) / Fs[s];
671 : d_cy_dz = d_faxial_dz[i + chainBeads * s] * (cos(2.0 * M_PI * Position[1]) - cy[s]) / Fs[s];
672 : d_ws_dz = (1 - pow(ws[s], 2)) * d_faxial_dz[i + chainBeads * s];
673 : d_Xsc_dz += (ws[s] * d_sx_dz + d_ws_dz * (sx[s] - Xsc)) / W;
674 : d_Xcc_dz += (ws[s] * d_cx_dz + d_ws_dz * (cx[s] - Xcc)) / W;
675 : d_Ysc_dz += (ws[s] * d_sy_dz + d_ws_dz * (sy[s] - Ysc)) / W;
676 : d_Ycc_dz += (ws[s] * d_cy_dz + d_ws_dz * (cy[s] - Ycc)) / W;
677 : }
678 : }
679 : d_Xcyl_dx[i] = auxX * (-Xsc * d_Xcc_dx + Xcc * d_Xsc_dx);
680 : d_Xcyl_dz[i] = auxX * (-Xsc * d_Xcc_dz + Xcc * d_Xsc_dz);
681 : d_Ycyl_dy[i] = auxY * (-Ysc * d_Ycc_dy + Ycc * d_Ysc_dy);
682 : d_Ycyl_dz[i] = auxY * (-Ysc * d_Ycc_dz + Ycc * d_Ysc_dz);
683 :
684 : ri = sqrt(pow(CylDistances[i][0], 2) + pow(CylDistances[i][1], 2));
685 : x = ri / RCYL[0];
686 : if (!((x <= -1.0 - HCH[0]) || (x >= 1.0 + HCH[0]))) {
687 : if (((-1.0 + HCH[0]) <= x) && (x <= (1.0 - HCH[0]))) {
688 : fradial = 1.0;
689 : } else if (((1.0 - HCH[0]) < x) && (x < (1.0 + HCH[0]))) {
690 : fradial = 0.5 - ((3.0 * x - 3.0) / (4.0 * HCH[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HCH[0], 3)));
691 : d_fradial_dx[i] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][0] / (RCYL[0] * ri);
692 : d_fradial_dy[i] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][1] / (RCYL[0] * ri);
693 : } else if (((-1.0 - HCH[0]) < x) && (x < (-1.0 + HCH[0]))) {
694 : fradial = 0.5 + ((3.0 * x + 3.0) / (4.0 * HCH[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HCH[0], 3)));
695 : d_fradial_dx[i] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][0] / (RCYL[0] * ri);
696 : d_fradial_dy[i] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][1] / (RCYL[0] * ri);
697 : }
698 :
699 : for (unsigned s = s1[i]; s <= s2[i]; s++) {
700 : Nsp[s] += fradial * faxial[i + chainBeads * s];
701 : Axs[s] += faxial[i + chainBeads * s] * d_fradial_dx[i];
702 : Ays[s] += faxial[i + chainBeads * s] * d_fradial_dy[i];
703 : fradial_d_faxial_dz[i + chainBeads * s] = fradial * d_faxial_dz[i + chainBeads * s];
704 : }
705 : }
706 : }
707 : }
708 :
709 184 : for (unsigned s = 0; s < NS[0]; s++) {
710 180 : if (Nsp[s] <= 1.0) {
711 149 : psi[s] = ZETA[0] * Nsp[s];
712 149 : d_psi[s] = ZETA[0];
713 149 : Xi_n += psi[s];
714 : } else {
715 31 : psi[s] = 1.0 - c * exp(-b * Nsp[s]);
716 31 : d_psi[s] = b * c * exp(-b * Nsp[s]);
717 31 : Xi_n += psi[s];
718 : }
719 : }
720 :
721 4 : Xi_n = Xi_n / NS[0];
722 :
723 : // Eq. 18 Hub & Awasthi JCTC 2017.
724 4 : std::vector<double> faxial_d_fradial_dx(chainBeads * NS[0]), faxial_d_fradial_dy(chainBeads * NS[0]), faxial_d_fradial_dz(chainBeads * NS[0]);
725 :
726 : // Eq. 13 Hub & Awasthi JCTC 2017 modified to considere the Heaviside_Chain step function (this only affect during the transition).
727 4 : std::vector<Vector> derivatives_Chain(chainBeads);
728 :
729 : #ifdef _OPENMP
730 : #if _OPENMP >= 201307
731 4 : #pragma omp parallel for private(aux)
732 : #endif
733 : #endif
734 : for (unsigned i = 0; i < chainBeads; i++) {
735 : if (analyzeThisParticle[i]) {
736 : for (unsigned s = s1[i]; s <= s2[i]; s++) {
737 : if (faxial[i + chainBeads * s]) {
738 : faxial_d_fradial_dx[i + chainBeads * s] = faxial[i + chainBeads * s] * d_fradial_dx[i] - d_Xcyl_dx[i] * Axs[s];
739 : faxial_d_fradial_dy[i + chainBeads * s] = faxial[i + chainBeads * s] * d_fradial_dy[i] - d_Ycyl_dy[i] * Ays[s];
740 : faxial_d_fradial_dz[i + chainBeads * s] = -d_Xcyl_dz[i] * Axs[s] - d_Ycyl_dz[i] * Ays[s];
741 : }
742 : }
743 :
744 : for (unsigned s = s1[i]; s <= s2[i]; s++) {
745 : aux = d_psi[s] / NS[0];
746 : derivatives_Chain[i][0] += aux * faxial_d_fradial_dx[i + chainBeads * s];
747 : derivatives_Chain[i][1] += aux * faxial_d_fradial_dy[i + chainBeads * s];
748 : derivatives_Chain[i][2] += aux * (faxial_d_fradial_dz[i + chainBeads * s] + fradial_d_faxial_dz[i + chainBeads * s]);
749 : }
750 : }
751 : }
752 :
753 : Tensor virial;
754 134604 : for (unsigned i = 0; i < chainBeads; i++) {
755 134600 : setAtomsDerivatives((i + noChainBeads), derivatives_Chain[i]);
756 269200 : virial -= Tensor(CylDistances[i], derivatives_Chain[i]);
757 : }
758 :
759 4 : setValue(Xi_n);
760 4 : setBoxDerivatives(virial);
761 4 : }
762 : }
763 : }
|
