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