File: | programs/mdrun/membed.c |
Location: | line 1195, column 9 |
Description: | Value stored to 'mem_nat' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by |
5 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
6 | * and including many others, as listed in the AUTHORS file in the |
7 | * top-level source directory and at http://www.gromacs.org. |
8 | * |
9 | * GROMACS is free software; you can redistribute it and/or |
10 | * modify it under the terms of the GNU Lesser General Public License |
11 | * as published by the Free Software Foundation; either version 2.1 |
12 | * of the License, or (at your option) any later version. |
13 | * |
14 | * GROMACS is distributed in the hope that it will be useful, |
15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 | * Lesser General Public License for more details. |
18 | * |
19 | * You should have received a copy of the GNU Lesser General Public |
20 | * License along with GROMACS; if not, see |
21 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
22 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
23 | * |
24 | * If you want to redistribute modifications to GROMACS, please |
25 | * consider that scientific software is very special. Version |
26 | * control is crucial - bugs must be traceable. We will be happy to |
27 | * consider code for inclusion in the official distribution, but |
28 | * derived work must not be called official GROMACS. Details are found |
29 | * in the README & COPYING files - if they are missing, get the |
30 | * official version at http://www.gromacs.org. |
31 | * |
32 | * To help us fund GROMACS development, we humbly ask that you cite |
33 | * the research papers on the package. Check out http://www.gromacs.org. |
34 | */ |
35 | #ifdef HAVE_CONFIG_H1 |
36 | #include <config.h> |
37 | #endif |
38 | |
39 | #include <signal.h> |
40 | #include <stdlib.h> |
41 | #include "typedefs.h" |
42 | #include "types/commrec.h" |
43 | #include "gromacs/utility/smalloc.h" |
44 | #include "gromacs/math/vec.h" |
45 | #include "macros.h" |
46 | #include "gromacs/utility/futil.h" |
47 | #include "gromacs/essentialdynamics/edsam.h" |
48 | #include "index.h" |
49 | #include "physics.h" |
50 | #include "names.h" |
51 | #include "mtop_util.h" |
52 | #include "gromacs/fileio/tpxio.h" |
53 | #include "gromacs/utility/cstringutil.h" |
54 | #include "membed.h" |
55 | #include "pbc.h" |
56 | #include "readinp.h" |
57 | #include "gromacs/gmxpreprocess/readir.h" |
58 | |
59 | /* information about scaling center */ |
60 | typedef struct { |
61 | rvec xmin; /* smallest coordinates of all embedded molecules */ |
62 | rvec xmax; /* largest coordinates of all embedded molecules */ |
63 | rvec *geom_cent; /* scaling center of each independent molecule to embed */ |
64 | int pieces; /* number of molecules to embed independently */ |
65 | int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */ |
66 | atom_id **subindex; /* atomids for independent molecule * |
67 | * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */ |
68 | } pos_ins_t; |
69 | |
70 | /* variables needed in do_md */ |
71 | struct membed { |
72 | int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */ |
73 | int it_z; /* same, but for z */ |
74 | real xy_step; /* stepsize used during resize in xy-plane */ |
75 | real z_step; /* same, but in z */ |
76 | rvec fac; /* initial scaling of the molecule to grow into the membrane */ |
77 | rvec *r_ins; /* final coordinates of the molecule to grow */ |
78 | pos_ins_t *pos_ins; /* scaling center for each piece to embed */ |
79 | }; |
80 | |
81 | /* membrane related variables */ |
82 | typedef struct { |
83 | char *name; /* name of index group to embed molecule into (usually membrane) */ |
84 | t_block mem_at; /* list all atoms in membrane */ |
85 | int nmol; /* number of membrane molecules overlapping with the molecule to embed */ |
86 | int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */ |
87 | real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/ |
88 | real zmin; /* minimum z coordinate of membrane */ |
89 | real zmax; /* maximum z coordinate of membrane */ |
90 | real zmed; /* median z coordinate of membrane */ |
91 | } mem_t; |
92 | |
93 | /* Lists all molecules in the membrane that overlap with the molecule to be embedded. * |
94 | * These will then be removed from the system */ |
95 | typedef struct { |
96 | int nr; /* number of molecules to remove */ |
97 | int *mol; /* list of molecule ids to remove */ |
98 | int *block; /* id of the molblock that the molecule to remove is part of */ |
99 | } rm_t; |
100 | |
101 | /* Get the global molecule id, and the corresponding molecule type and id of the * |
102 | * molblock from the global atom nr. */ |
103 | static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block) |
104 | { |
105 | int mol_id = 0; |
106 | int i; |
107 | int atnr_mol; |
108 | gmx_mtop_atomlookup_t alook; |
109 | |
110 | alook = gmx_mtop_atomlookup_settle_init(mtop); |
111 | gmx_mtop_atomnr_to_molblock_ind(alook, at, block, &mol_id, &atnr_mol); |
112 | for (i = 0; i < *block; i++) |
113 | { |
114 | mol_id += mtop->molblock[i].nmol; |
115 | } |
116 | *type = mtop->molblock[*block].type; |
117 | |
118 | gmx_mtop_atomlookup_destroy(alook); |
119 | |
120 | return mol_id; |
121 | } |
122 | |
123 | /* Get the id of the molblock from a global molecule id */ |
124 | static int get_molblock(int mol_id, int nmblock, gmx_molblock_t *mblock) |
125 | { |
126 | int i; |
127 | int nmol = 0; |
128 | |
129 | for (i = 0; i < nmblock; i++) |
130 | { |
131 | nmol += mblock[i].nmol; |
132 | if (mol_id < nmol) |
133 | { |
134 | return i; |
135 | } |
136 | } |
137 | |
138 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 138, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol); |
139 | |
140 | return -1; |
141 | } |
142 | |
143 | static int get_tpr_version(const char *infile) |
144 | { |
145 | t_tpxheader header; |
146 | int version, generation; |
147 | |
148 | read_tpxheader(infile, &header, TRUE1, &version, &generation); |
149 | |
150 | return version; |
151 | } |
152 | |
153 | /* Get a list of all the molecule types that are present in a group of atoms. * |
154 | * Because all interaction within the group to embed are removed on the topology * |
155 | * level, if the same molecule type is found in another part of the system, these * |
156 | * would also be affected. Therefore we have to check if the embedded and rest group * |
157 | * share common molecule types. If so, membed will stop with an error. */ |
158 | static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist) |
159 | { |
160 | int i, j, nr, mol_id; |
161 | int type = 0, block = 0; |
162 | gmx_bool bNEW; |
163 | |
164 | nr = 0; |
165 | snew(tlist->index, at->nr)(tlist->index) = save_calloc("tlist->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 165, (at->nr), sizeof(*(tlist->index))); |
166 | for (i = 0; i < at->nr; i++) |
167 | { |
168 | bNEW = TRUE1; |
169 | mol_id = get_mol_id(at->index[i], mtop, &type, &block); |
170 | for (j = 0; j < nr; j++) |
171 | { |
172 | if (tlist->index[j] == type) |
173 | { |
174 | bNEW = FALSE0; |
175 | } |
176 | } |
177 | |
178 | if (bNEW == TRUE1) |
179 | { |
180 | tlist->index[nr] = type; |
181 | nr++; |
182 | } |
183 | } |
184 | srenew(tlist->index, nr)(tlist->index) = save_realloc("tlist->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 184, (tlist->index), (nr), sizeof(*(tlist->index))); |
185 | return nr; |
186 | } |
187 | |
188 | /* Do the actual check of the molecule types between embedded and rest group */ |
189 | static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop) |
190 | { |
191 | t_block *ins_mtype, *rest_mtype; |
192 | int i, j; |
193 | |
194 | snew(ins_mtype, 1)(ins_mtype) = save_calloc("ins_mtype", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 194, (1), sizeof(*(ins_mtype))); |
195 | snew(rest_mtype, 1)(rest_mtype) = save_calloc("rest_mtype", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 195, (1), sizeof(*(rest_mtype))); |
196 | ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype ); |
197 | rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype); |
198 | |
199 | for (i = 0; i < ins_mtype->nr; i++) |
200 | { |
201 | for (j = 0; j < rest_mtype->nr; j++) |
202 | { |
203 | if (ins_mtype->index[i] == rest_mtype->index[j]) |
204 | { |
205 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 205, "Moleculetype %s is found both in the group to insert and the rest of the system.\n" |
206 | "1. Your *.ndx and *.top do not match\n" |
207 | "2. You are inserting some molecules of type %s (for example xray-solvent), while\n" |
208 | "the same moleculetype is also used in the rest of the system (solvent box). Because\n" |
209 | "we need to exclude all interactions between the atoms in the group to\n" |
210 | "insert, the same moleculetype can not be used in both groups. Change the\n" |
211 | "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n" |
212 | "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name), |
213 | *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name)); |
214 | } |
215 | } |
216 | } |
217 | |
218 | sfree(ins_mtype->index)save_free("ins_mtype->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 218, (ins_mtype->index)); |
219 | sfree(rest_mtype->index)save_free("rest_mtype->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 219, (rest_mtype->index)); |
220 | sfree(ins_mtype)save_free("ins_mtype", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 220, (ins_mtype)); |
221 | sfree(rest_mtype)save_free("rest_mtype", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 221, (rest_mtype)); |
222 | } |
223 | |
224 | static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max, |
225 | int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn, |
226 | int *pieces, gmx_bool *bALLOW_ASYMMETRY) |
227 | { |
228 | warninp_t wi; |
229 | t_inpfile *inp; |
230 | int ninp; |
231 | |
232 | wi = init_warning(TRUE1, 0); |
233 | |
234 | inp = read_inpfile(membed_input, &ninp, wi); |
235 | ITYPE ("nxy", *it_xy, 1000)*it_xy = get_eint(&ninp, &inp, "nxy", 1000, wi); |
236 | ITYPE ("nz", *it_z, 0)*it_z = get_eint(&ninp, &inp, "nz", 0, wi); |
237 | RTYPE ("xyinit", *xy_fac, 0.5)*xy_fac = get_ereal(&ninp, &inp, "xyinit", 0.5, wi); |
238 | RTYPE ("xyend", *xy_max, 1.0)*xy_max = get_ereal(&ninp, &inp, "xyend", 1.0, wi); |
239 | RTYPE ("zinit", *z_fac, 1.0)*z_fac = get_ereal(&ninp, &inp, "zinit", 1.0, wi); |
240 | RTYPE ("zend", *z_max, 1.0)*z_max = get_ereal(&ninp, &inp, "zend", 1.0, wi); |
241 | RTYPE ("rad", *probe_rad, 0.22)*probe_rad = get_ereal(&ninp, &inp, "rad", 0.22, wi); |
242 | ITYPE ("ndiff", *low_up_rm, 0)*low_up_rm = get_eint(&ninp, &inp, "ndiff", 0, wi); |
243 | ITYPE ("maxwarn", *maxwarn, 0)*maxwarn = get_eint(&ninp, &inp, "maxwarn", 0, wi); |
244 | ITYPE ("pieces", *pieces, 1)*pieces = get_eint(&ninp, &inp, "pieces", 1, wi); |
245 | EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names)*bALLOW_ASYMMETRY = get_eeenum(&ninp, &inp, "asymmetry" , yesno_names, wi); |
246 | |
247 | write_inpfile(membed_input, ninp, inp, FALSE0, wi); |
248 | } |
249 | |
250 | /* Obtain the maximum and minimum coordinates of the group to be embedded */ |
251 | static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins, |
252 | gmx_groups_t *groups, int ins_grp_id, real xy_max) |
253 | { |
254 | int i, gid, c = 0; |
255 | real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax; |
256 | const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to * |
257 | * determine the overlap between molecule to embed and membrane */ |
258 | const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */ |
259 | snew(rest_at->index, state->natoms)(rest_at->index) = save_calloc("rest_at->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 259, (state->natoms), sizeof(*(rest_at->index))); |
260 | |
261 | xmin = xmax = state->x[ins_at->index[0]][XX0]; |
262 | ymin = ymax = state->x[ins_at->index[0]][YY1]; |
263 | zmin = zmax = state->x[ins_at->index[0]][ZZ2]; |
264 | |
265 | for (i = 0; i < state->natoms; i++) |
266 | { |
267 | gid = groups->grpnr[egcFREEZE][i]; |
268 | if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id) |
269 | { |
270 | x = state->x[i][XX0]; |
271 | if (x < xmin) |
272 | { |
273 | xmin = x; |
274 | } |
275 | if (x > xmax) |
276 | { |
277 | xmax = x; |
278 | } |
279 | y = state->x[i][YY1]; |
280 | if (y < ymin) |
281 | { |
282 | ymin = y; |
283 | } |
284 | if (y > ymax) |
285 | { |
286 | ymax = y; |
287 | } |
288 | z = state->x[i][ZZ2]; |
289 | if (z < zmin) |
290 | { |
291 | zmin = z; |
292 | } |
293 | if (z > zmax) |
294 | { |
295 | zmax = z; |
296 | } |
297 | } |
298 | else |
299 | { |
300 | rest_at->index[c] = i; |
301 | c++; |
302 | } |
303 | } |
304 | |
305 | rest_at->nr = c; |
306 | srenew(rest_at->index, c)(rest_at->index) = save_realloc("rest_at->index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 306, (rest_at->index), (c), sizeof(*(rest_at->index)) ); |
307 | |
308 | if (xy_max > fac_inp_size) |
309 | { |
310 | pos_ins->xmin[XX0] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2; |
311 | pos_ins->xmin[YY1] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2; |
312 | |
313 | pos_ins->xmax[XX0] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2; |
314 | pos_ins->xmax[YY1] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2; |
315 | } |
316 | else |
317 | { |
318 | pos_ins->xmin[XX0] = xmin; |
319 | pos_ins->xmin[YY1] = ymin; |
320 | |
321 | pos_ins->xmax[XX0] = xmax; |
322 | pos_ins->xmax[YY1] = ymax; |
323 | } |
324 | |
325 | if ( (zmax-zmin) < min_memthick) |
326 | { |
327 | pos_ins->xmin[ZZ2] = zmin+(zmax-zmin)/2.0-0.5*min_memthick; |
328 | pos_ins->xmax[ZZ2] = zmin+(zmax-zmin)/2.0+0.5*min_memthick; |
329 | } |
330 | else |
331 | { |
332 | pos_ins->xmin[ZZ2] = zmin; |
333 | pos_ins->xmax[ZZ2] = zmax; |
334 | } |
335 | |
336 | return c; |
337 | } |
338 | |
339 | /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the * |
340 | * xy-plane and counting the number of occupied grid points */ |
341 | static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p) |
342 | { |
343 | real x, y, dx = 0.15, dy = 0.15, area = 0.0; |
344 | real add, memmin, memmax; |
345 | int c, at; |
346 | |
347 | /* min and max membrane coordinate are altered to reduce the influence of the * |
348 | * boundary region */ |
349 | memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin); |
350 | memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin); |
351 | |
352 | for (x = pos_ins->xmin[XX0]; x < pos_ins->xmax[XX0]; x += dx) |
353 | { |
354 | for (y = pos_ins->xmin[YY1]; y < pos_ins->xmax[YY1]; y += dy) |
355 | { |
356 | c = 0; |
357 | add = 0.0; |
358 | do |
359 | { |
360 | at = ins_at->index[c]; |
361 | if ( (r[at][XX0] >= x) && (r[at][XX0] < x+dx) && |
362 | (r[at][YY1] >= y) && (r[at][YY1] < y+dy) && |
363 | (r[at][ZZ2] > memmin) && (r[at][ZZ2] < memmax) ) |
364 | { |
365 | add = 1.0; |
366 | } |
367 | c++; |
368 | } |
369 | while ( (c < ins_at->nr) && (add < 0.5) ); |
370 | area += add; |
371 | } |
372 | } |
373 | area = area*dx*dy; |
374 | |
375 | return area; |
376 | } |
377 | |
378 | static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins) |
379 | { |
380 | int i, j, at, mol, nmol, nmolbox, count; |
381 | t_block *mem_a; |
382 | real z, zmin, zmax, mem_area; |
383 | gmx_bool bNew; |
384 | atom_id *mol_id; |
385 | int type = 0, block = 0; |
386 | |
387 | nmol = count = 0; |
388 | mem_a = &(mem_p->mem_at); |
389 | snew(mol_id, mem_a->nr)(mol_id) = save_calloc("mol_id", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 389, (mem_a->nr), sizeof(*(mol_id))); |
390 | zmin = pos_ins->xmax[ZZ2]; |
391 | zmax = pos_ins->xmin[ZZ2]; |
392 | for (i = 0; i < mem_a->nr; i++) |
393 | { |
394 | at = mem_a->index[i]; |
395 | if ( (r[at][XX0] > pos_ins->xmin[XX0]) && (r[at][XX0] < pos_ins->xmax[XX0]) && |
396 | (r[at][YY1] > pos_ins->xmin[YY1]) && (r[at][YY1] < pos_ins->xmax[YY1]) && |
397 | (r[at][ZZ2] > pos_ins->xmin[ZZ2]) && (r[at][ZZ2] < pos_ins->xmax[ZZ2]) ) |
398 | { |
399 | mol = get_mol_id(at, mtop, &type, &block); |
400 | bNew = TRUE1; |
401 | for (j = 0; j < nmol; j++) |
402 | { |
403 | if (mol == mol_id[j]) |
404 | { |
405 | bNew = FALSE0; |
406 | } |
407 | } |
408 | |
409 | if (bNew) |
410 | { |
411 | mol_id[nmol] = mol; |
412 | nmol++; |
413 | } |
414 | |
415 | z = r[at][ZZ2]; |
416 | if (z < zmin) |
417 | { |
418 | zmin = z; |
419 | } |
420 | |
421 | if (z > zmax) |
422 | { |
423 | zmax = z; |
424 | } |
425 | |
426 | count++; |
427 | } |
428 | } |
429 | |
430 | mem_p->nmol = nmol; |
431 | srenew(mol_id, nmol)(mol_id) = save_realloc("mol_id", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 431, (mol_id), (nmol), sizeof(*(mol_id))); |
432 | mem_p->mol_id = mol_id; |
433 | |
434 | if ((zmax-zmin) > (box[ZZ2][ZZ2]-0.5)) |
435 | { |
436 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 436, "Something is wrong with your membrane. Max and min z values are %f and %f.\n" |
437 | "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n" |
438 | "so that one membrane is distributed over two periodic box images. Another possibility is that\n" |
439 | "your water layer is not thick enough.\n", zmax, zmin); |
440 | } |
441 | mem_p->zmin = zmin; |
442 | mem_p->zmax = zmax; |
443 | mem_p->zmed = (zmax-zmin)/2+zmin; |
444 | |
445 | /*number of membrane molecules in protein box*/ |
446 | nmolbox = count/mtop->molblock[block].natoms_mol; |
447 | /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/ |
448 | mem_area = (pos_ins->xmax[XX0]-pos_ins->xmin[XX0])*(pos_ins->xmax[YY1]-pos_ins->xmin[YY1]); |
449 | /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/ |
450 | mem_p->lip_area = 2.0*mem_area/(double)nmolbox; |
451 | |
452 | return mem_p->mem_at.nr; |
453 | } |
454 | |
455 | static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r, |
456 | gmx_bool bALLOW_ASYMMETRY) |
457 | { |
458 | int i, j, at, c, outsidesum, gctr = 0; |
459 | int idxsum = 0; |
460 | |
461 | /*sanity check*/ |
462 | for (i = 0; i < pos_ins->pieces; i++) |
463 | { |
464 | idxsum += pos_ins->nidx[i]; |
465 | } |
466 | |
467 | if (idxsum != ins_at->nr) |
468 | { |
469 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 469, "Piecewise sum of inserted atoms not same as size of group selected to insert."); |
470 | } |
471 | |
472 | snew(pos_ins->geom_cent, pos_ins->pieces)(pos_ins->geom_cent) = save_calloc("pos_ins->geom_cent" , "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c", 472, (pos_ins->pieces), sizeof(*(pos_ins->geom_cent))); |
473 | for (i = 0; i < pos_ins->pieces; i++) |
474 | { |
475 | c = 0; |
476 | outsidesum = 0; |
477 | for (j = 0; j < DIM3; j++) |
478 | { |
479 | pos_ins->geom_cent[i][j] = 0; |
480 | } |
481 | |
482 | for (j = 0; j < pos_ins->nidx[i]; j++) |
483 | { |
484 | at = pos_ins->subindex[i][j]; |
485 | copy_rvec(r[at], r_ins[gctr]); |
486 | if ( (r_ins[gctr][ZZ2] < mem_p->zmax) && (r_ins[gctr][ZZ2] > mem_p->zmin) ) |
487 | { |
488 | rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]); |
489 | c++; |
490 | } |
491 | else |
492 | { |
493 | outsidesum++; |
494 | } |
495 | gctr++; |
496 | } |
497 | |
498 | if (c > 0) |
499 | { |
500 | svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]); |
501 | } |
502 | |
503 | if (!bALLOW_ASYMMETRY) |
504 | { |
505 | pos_ins->geom_cent[i][ZZ2] = mem_p->zmed; |
506 | } |
507 | |
508 | fprintf(stderrstderr, "Embedding piece %d with center of geometry: %f %f %f\n", |
509 | i, pos_ins->geom_cent[i][XX0], pos_ins->geom_cent[i][YY1], pos_ins->geom_cent[i][ZZ2]); |
510 | } |
511 | fprintf(stderrstderr, "\n"); |
512 | } |
513 | |
514 | /* resize performed in the md loop */ |
515 | static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, rvec fac) |
516 | { |
517 | int i, j, k, at, c = 0; |
518 | for (k = 0; k < pos_ins->pieces; k++) |
519 | { |
520 | for (i = 0; i < pos_ins->nidx[k]; i++) |
521 | { |
522 | at = pos_ins->subindex[k][i]; |
523 | for (j = 0; j < DIM3; j++) |
524 | { |
525 | r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]); |
526 | } |
527 | c++; |
528 | } |
529 | } |
530 | } |
531 | |
532 | /* generate the list of membrane molecules that overlap with the molecule to be embedded. * |
533 | * The molecule to be embedded is already reduced in size. */ |
534 | static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop, |
535 | rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, |
536 | int low_up_rm, gmx_bool bALLOW_ASYMMETRY) |
537 | { |
538 | int i, j, k, l, at, at2, mol_id; |
539 | int type = 0, block = 0; |
540 | int nrm, nupper, nlower; |
541 | real r_min_rad, z_lip, min_norm; |
542 | gmx_bool bRM; |
543 | rvec dr, dr_tmp; |
544 | real *dist; |
545 | int *order; |
546 | |
547 | r_min_rad = probe_rad*probe_rad; |
548 | snew(rm_p->mol, mtop->mols.nr)(rm_p->mol) = save_calloc("rm_p->mol", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 548, (mtop->mols.nr), sizeof(*(rm_p->mol))); |
549 | snew(rm_p->block, mtop->mols.nr)(rm_p->block) = save_calloc("rm_p->block", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 549, (mtop->mols.nr), sizeof(*(rm_p->block))); |
550 | nrm = nupper = 0; |
551 | nlower = low_up_rm; |
552 | for (i = 0; i < ins_at->nr; i++) |
553 | { |
554 | at = ins_at->index[i]; |
555 | for (j = 0; j < rest_at->nr; j++) |
556 | { |
557 | at2 = rest_at->index[j]; |
558 | pbc_dx(pbc, r[at], r[at2], dr); |
559 | |
560 | if (norm2(dr) < r_min_rad) |
561 | { |
562 | mol_id = get_mol_id(at2, mtop, &type, &block); |
563 | bRM = TRUE1; |
564 | for (l = 0; l < nrm; l++) |
565 | { |
566 | if (rm_p->mol[l] == mol_id) |
567 | { |
568 | bRM = FALSE0; |
569 | } |
570 | } |
571 | |
572 | if (bRM) |
573 | { |
574 | rm_p->mol[nrm] = mol_id; |
575 | rm_p->block[nrm] = block; |
576 | nrm++; |
577 | z_lip = 0.0; |
578 | for (l = 0; l < mem_p->nmol; l++) |
579 | { |
580 | if (mol_id == mem_p->mol_id[l]) |
581 | { |
582 | for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++) |
583 | { |
584 | z_lip += r[k][ZZ2]; |
585 | } |
586 | z_lip /= mtop->molblock[block].natoms_mol; |
587 | if (z_lip < mem_p->zmed) |
588 | { |
589 | nlower++; |
590 | } |
591 | else |
592 | { |
593 | nupper++; |
594 | } |
595 | } |
596 | } |
597 | } |
598 | } |
599 | } |
600 | } |
601 | |
602 | /*make sure equal number of lipids from upper and lower layer are removed */ |
603 | if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) ) |
604 | { |
605 | snew(dist, mem_p->nmol)(dist) = save_calloc("dist", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 605, (mem_p->nmol), sizeof(*(dist))); |
606 | snew(order, mem_p->nmol)(order) = save_calloc("order", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 606, (mem_p->nmol), sizeof(*(order))); |
607 | for (i = 0; i < mem_p->nmol; i++) |
608 | { |
609 | at = mtop->mols.index[mem_p->mol_id[i]]; |
610 | pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr); |
611 | if (pos_ins->pieces > 1) |
612 | { |
613 | /*minimum dr value*/ |
614 | min_norm = norm2(dr); |
615 | for (k = 1; k < pos_ins->pieces; k++) |
616 | { |
617 | pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp); |
618 | if (norm2(dr_tmp) < min_norm) |
619 | { |
620 | min_norm = norm2(dr_tmp); |
621 | copy_rvec(dr_tmp, dr); |
622 | } |
623 | } |
624 | } |
625 | dist[i] = dr[XX0]*dr[XX0]+dr[YY1]*dr[YY1]; |
626 | j = i-1; |
627 | while (j >= 0 && dist[i] < dist[order[j]]) |
628 | { |
629 | order[j+1] = order[j]; |
630 | j--; |
631 | } |
632 | order[j+1] = i; |
633 | } |
634 | |
635 | i = 0; |
636 | while (nupper != nlower) |
637 | { |
638 | mol_id = mem_p->mol_id[order[i]]; |
639 | block = get_molblock(mol_id, mtop->nmolblock, mtop->molblock); |
640 | bRM = TRUE1; |
641 | for (l = 0; l < nrm; l++) |
642 | { |
643 | if (rm_p->mol[l] == mol_id) |
644 | { |
645 | bRM = FALSE0; |
646 | } |
647 | } |
648 | |
649 | if (bRM) |
650 | { |
651 | z_lip = 0; |
652 | for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++) |
653 | { |
654 | z_lip += r[k][ZZ2]; |
655 | } |
656 | z_lip /= mtop->molblock[block].natoms_mol; |
657 | if (nupper > nlower && z_lip < mem_p->zmed) |
658 | { |
659 | rm_p->mol[nrm] = mol_id; |
660 | rm_p->block[nrm] = block; |
661 | nrm++; |
662 | nlower++; |
663 | } |
664 | else if (nupper < nlower && z_lip > mem_p->zmed) |
665 | { |
666 | rm_p->mol[nrm] = mol_id; |
667 | rm_p->block[nrm] = block; |
668 | nrm++; |
669 | nupper++; |
670 | } |
671 | } |
672 | i++; |
673 | |
674 | if (i > mem_p->nmol) |
675 | { |
676 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 676, "Trying to remove more lipid molecules than there are in the membrane"); |
677 | } |
678 | } |
679 | sfree(dist)save_free("dist", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 679, (dist)); |
680 | sfree(order)save_free("order", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 680, (order)); |
681 | } |
682 | |
683 | rm_p->nr = nrm; |
684 | srenew(rm_p->mol, nrm)(rm_p->mol) = save_realloc("rm_p->mol", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 684, (rm_p->mol), (nrm), sizeof(*(rm_p->mol))); |
685 | srenew(rm_p->block, nrm)(rm_p->block) = save_realloc("rm_p->block", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 685, (rm_p->block), (nrm), sizeof(*(rm_p->block))); |
686 | |
687 | return nupper+nlower; |
688 | } |
689 | |
690 | /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/ |
691 | static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state, |
692 | t_block *ins_at, pos_ins_t *pos_ins) |
693 | { |
694 | int i, j, k, n, rm, mol_id, at, block; |
695 | rvec *x_tmp, *v_tmp; |
696 | atom_id *list, *new_mols; |
697 | unsigned char *new_egrp[egcNR]; |
698 | gmx_bool bRM; |
699 | int RMmolblock; |
700 | |
701 | snew(list, state->natoms)(list) = save_calloc("list", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 701, (state->natoms), sizeof(*(list))); |
702 | n = 0; |
703 | for (i = 0; i < rm_p->nr; i++) |
704 | { |
705 | mol_id = rm_p->mol[i]; |
706 | at = mtop->mols.index[mol_id]; |
707 | block = rm_p->block[i]; |
708 | mtop->molblock[block].nmol--; |
709 | for (j = 0; j < mtop->molblock[block].natoms_mol; j++) |
710 | { |
711 | list[n] = at+j; |
712 | n++; |
713 | } |
714 | mtop->mols.index[mol_id] = -1; |
715 | } |
716 | |
717 | mtop->mols.nr -= rm_p->nr; |
718 | mtop->mols.nalloc_index -= rm_p->nr; |
719 | snew(new_mols, mtop->mols.nr)(new_mols) = save_calloc("new_mols", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 719, (mtop->mols.nr), sizeof(*(new_mols))); |
720 | for (i = 0; i < mtop->mols.nr+rm_p->nr; i++) |
721 | { |
722 | j = 0; |
723 | if (mtop->mols.index[i] != -1) |
724 | { |
725 | new_mols[j] = mtop->mols.index[i]; |
726 | j++; |
727 | } |
728 | } |
729 | sfree(mtop->mols.index)save_free("mtop->mols.index", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 729, (mtop->mols.index)); |
730 | mtop->mols.index = new_mols; |
731 | mtop->natoms -= n; |
732 | state->natoms -= n; |
733 | state->nalloc = state->natoms; |
734 | snew(x_tmp, state->nalloc)(x_tmp) = save_calloc("x_tmp", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 734, (state->nalloc), sizeof(*(x_tmp))); |
735 | snew(v_tmp, state->nalloc)(v_tmp) = save_calloc("v_tmp", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 735, (state->nalloc), sizeof(*(v_tmp))); |
736 | |
737 | for (i = 0; i < egcNR; i++) |
738 | { |
739 | if (groups->grpnr[i] != NULL((void*)0)) |
740 | { |
741 | groups->ngrpnr[i] = state->natoms; |
742 | snew(new_egrp[i], state->natoms)(new_egrp[i]) = save_calloc("new_egrp[i]", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 742, (state->natoms), sizeof(*(new_egrp[i]))); |
743 | } |
744 | } |
745 | |
746 | rm = 0; |
747 | for (i = 0; i < state->natoms+n; i++) |
748 | { |
749 | bRM = FALSE0; |
750 | for (j = 0; j < n; j++) |
751 | { |
752 | if (i == list[j]) |
753 | { |
754 | bRM = TRUE1; |
755 | rm++; |
756 | } |
757 | } |
758 | |
759 | if (!bRM) |
760 | { |
761 | for (j = 0; j < egcNR; j++) |
762 | { |
763 | if (groups->grpnr[j] != NULL((void*)0)) |
764 | { |
765 | new_egrp[j][i-rm] = groups->grpnr[j][i]; |
766 | } |
767 | } |
768 | copy_rvec(state->x[i], x_tmp[i-rm]); |
769 | copy_rvec(state->v[i], v_tmp[i-rm]); |
770 | for (j = 0; j < ins_at->nr; j++) |
771 | { |
772 | if (i == ins_at->index[j]) |
773 | { |
774 | ins_at->index[j] = i-rm; |
775 | } |
776 | } |
777 | |
778 | for (j = 0; j < pos_ins->pieces; j++) |
779 | { |
780 | for (k = 0; k < pos_ins->nidx[j]; k++) |
781 | { |
782 | if (i == pos_ins->subindex[j][k]) |
783 | { |
784 | pos_ins->subindex[j][k] = i-rm; |
785 | } |
786 | } |
787 | } |
788 | } |
789 | } |
790 | sfree(state->x)save_free("state->x", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 790, (state->x)); |
791 | state->x = x_tmp; |
792 | sfree(state->v)save_free("state->v", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 792, (state->v)); |
793 | state->v = v_tmp; |
794 | |
795 | for (i = 0; i < egcNR; i++) |
796 | { |
797 | if (groups->grpnr[i] != NULL((void*)0)) |
798 | { |
799 | sfree(groups->grpnr[i])save_free("groups->grpnr[i]", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 799, (groups->grpnr[i])); |
800 | groups->grpnr[i] = new_egrp[i]; |
801 | } |
802 | } |
803 | |
804 | /* remove empty molblocks */ |
805 | RMmolblock = 0; |
806 | for (i = 0; i < mtop->nmolblock; i++) |
807 | { |
808 | if (mtop->molblock[i].nmol == 0) |
809 | { |
810 | RMmolblock++; |
811 | } |
812 | else |
813 | { |
814 | mtop->molblock[i-RMmolblock] = mtop->molblock[i]; |
815 | } |
816 | } |
817 | mtop->nmolblock -= RMmolblock; |
818 | } |
819 | |
820 | /* remove al bonded interactions from mtop for the molecule to be embedded */ |
821 | int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop) |
822 | { |
823 | int i, j, m; |
824 | int type, natom, nmol, at, atom1 = 0, rm_at = 0; |
825 | gmx_bool *bRM, bINS; |
826 | /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/ |
827 | /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make * |
828 | * sure that g_membed exits with a warning when there are molecules of the same type not in the * |
829 | * ins_at index group. MGWolf 050710 */ |
830 | |
831 | |
832 | snew(bRM, mtop->nmoltype)(bRM) = save_calloc("bRM", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 832, (mtop->nmoltype), sizeof(*(bRM))); |
833 | for (i = 0; i < mtop->nmoltype; i++) |
834 | { |
835 | bRM[i] = TRUE1; |
836 | } |
837 | |
838 | for (i = 0; i < mtop->nmolblock; i++) |
839 | { |
840 | /*loop over molecule blocks*/ |
841 | type = mtop->molblock[i].type; |
842 | natom = mtop->molblock[i].natoms_mol; |
843 | nmol = mtop->molblock[i].nmol; |
844 | |
845 | for (j = 0; j < natom*nmol && bRM[type] == TRUE1; j++) |
846 | { |
847 | /*loop over atoms in the block*/ |
848 | at = j+atom1; /*atom index = block index + offset*/ |
849 | bINS = FALSE0; |
850 | |
851 | for (m = 0; (m < ins_at->nr) && (bINS == FALSE0); m++) |
852 | { |
853 | /*loop over atoms in insertion index group to determine if we're inserting one*/ |
854 | if (at == ins_at->index[m]) |
855 | { |
856 | bINS = TRUE1; |
857 | } |
858 | } |
859 | bRM[type] = bINS; |
860 | } |
861 | atom1 += natom*nmol; /*update offset*/ |
862 | if (bRM[type]) |
863 | { |
864 | rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/ |
865 | } |
866 | } |
867 | |
868 | for (i = 0; i < mtop->nmoltype; i++) |
869 | { |
870 | if (bRM[i]) |
871 | { |
872 | for (j = 0; j < F_LJ; j++) |
873 | { |
874 | mtop->moltype[i].ilist[j].nr = 0; |
875 | } |
876 | |
877 | for (j = F_POSRES; j <= F_VSITEN; j++) |
878 | { |
879 | mtop->moltype[i].ilist[j].nr = 0; |
880 | } |
881 | } |
882 | } |
883 | sfree(bRM)save_free("bRM", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 883, (bRM)); |
884 | |
885 | return rm_at; |
886 | } |
887 | |
888 | /* Write a topology where the number of molecules is correct for the system after embedding */ |
889 | static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop) |
890 | { |
891 | #define TEMP_FILENM "temp.top" |
892 | int bMolecules = 0; |
893 | FILE *fpin, *fpout; |
894 | char buf[STRLEN4096], buf2[STRLEN4096], *temp; |
895 | int i, *nmol_rm, nmol, line; |
896 | |
897 | fpin = gmx_ffopen(topfile, "r"); |
898 | fpout = gmx_ffopen(TEMP_FILENM, "w"); |
899 | |
900 | snew(nmol_rm, mtop->nmoltype)(nmol_rm) = save_calloc("nmol_rm", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 900, (mtop->nmoltype), sizeof(*(nmol_rm))); |
901 | for (i = 0; i < rm_p->nr; i++) |
902 | { |
903 | nmol_rm[rm_p->block[i]]++; |
904 | } |
905 | |
906 | line = 0; |
907 | while (fgets(buf, STRLEN4096, fpin)) |
908 | { |
909 | line++; |
910 | if (buf[0] != ';') |
911 | { |
912 | strcpy(buf2, buf); |
913 | if ((temp = strchr(buf2, '\n')(__extension__ (__builtin_constant_p ('\n') && !__builtin_constant_p (buf2) && ('\n') == '\0' ? (char *) __rawmemchr (buf2 , '\n') : __builtin_strchr (buf2, '\n')))) != NULL((void*)0)) |
914 | { |
915 | temp[0] = '\0'; |
916 | } |
917 | ltrim(buf2); |
918 | if (buf2[0] == '[') |
919 | { |
920 | buf2[0] = ' '; |
921 | if ((temp = strchr(buf2, '\n')(__extension__ (__builtin_constant_p ('\n') && !__builtin_constant_p (buf2) && ('\n') == '\0' ? (char *) __rawmemchr (buf2 , '\n') : __builtin_strchr (buf2, '\n')))) != NULL((void*)0)) |
922 | { |
923 | temp[0] = '\0'; |
924 | } |
925 | rtrim(buf2); |
926 | if (buf2[strlen(buf2)-1] == ']') |
927 | { |
928 | buf2[strlen(buf2)-1] = '\0'; |
929 | ltrim(buf2); |
930 | rtrim(buf2); |
931 | if (gmx_strcasecmp(buf2, "molecules") == 0) |
932 | { |
933 | bMolecules = 1; |
934 | } |
935 | } |
936 | fprintf(fpout, "%s", buf); |
937 | } |
938 | else if (bMolecules == 1) |
939 | { |
940 | for (i = 0; i < mtop->nmolblock; i++) |
941 | { |
942 | nmol = mtop->molblock[i].nmol; |
943 | sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol); |
944 | fprintf(fpout, "%s", buf); |
945 | } |
946 | bMolecules = 2; |
947 | } |
948 | else if (bMolecules == 2) |
949 | { |
950 | /* print nothing */ |
951 | } |
952 | else |
953 | { |
954 | fprintf(fpout, "%s", buf); |
955 | } |
956 | } |
957 | else |
958 | { |
959 | fprintf(fpout, "%s", buf); |
960 | } |
961 | } |
962 | |
963 | gmx_ffclose(fpout); |
964 | /* use gmx_ffopen to generate backup of topinout */ |
965 | fpout = gmx_ffopen(topfile, "w"); |
966 | gmx_ffclose(fpout); |
967 | rename(TEMP_FILENM, topfile); |
968 | #undef TEMP_FILENM |
969 | } |
970 | |
971 | void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x) |
972 | { |
973 | /* Set new positions for the group to embed */ |
974 | if (step_rel <= membed->it_xy) |
975 | { |
976 | membed->fac[0] += membed->xy_step; |
977 | membed->fac[1] += membed->xy_step; |
978 | } |
979 | else if (step_rel <= (membed->it_xy+membed->it_z)) |
980 | { |
981 | membed->fac[2] += membed->z_step; |
982 | } |
983 | resize(membed->r_ins, x, membed->pos_ins, membed->fac); |
984 | } |
985 | |
986 | gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop, |
987 | t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt) |
988 | { |
989 | char *ins, **gnames; |
990 | int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0; |
991 | int ng, j, max_lip_rm, ins_grp_id, ins_nat, mem_nat, ntype, lip_rm, tpr_version; |
992 | real prot_area; |
993 | rvec *r_ins = NULL((void*)0); |
994 | t_block *ins_at, *rest_at; |
995 | pos_ins_t *pos_ins; |
996 | mem_t *mem_p; |
997 | rm_t *rm_p; |
998 | gmx_groups_t *groups; |
999 | gmx_bool bExcl = FALSE0; |
1000 | t_atoms atoms; |
1001 | t_pbc *pbc; |
1002 | char **piecename = NULL((void*)0); |
1003 | gmx_membed_t membed = NULL((void*)0); |
1004 | |
1005 | /* input variables */ |
1006 | const char *membed_input; |
1007 | real xy_fac = 0.5; |
1008 | real xy_max = 1.0; |
1009 | real z_fac = 1.0; |
1010 | real z_max = 1.0; |
1011 | int it_xy = 1000; |
1012 | int it_z = 0; |
1013 | real probe_rad = 0.22; |
1014 | int low_up_rm = 0; |
1015 | int maxwarn = 0; |
1016 | int pieces = 1; |
1017 | gmx_bool bALLOW_ASYMMETRY = FALSE0; |
1018 | |
1019 | /* sanity check constants */ /* Issue a warning when: */ |
1020 | const int membed_version = 58; /* tpr version is smaller */ |
1021 | const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule * |
1022 | * and rest smaller than this value is probably too small */ |
1023 | const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */ |
1024 | const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */ |
1025 | const int min_it_z = 100; /* the number of steps to embed in z is smaller */ |
1026 | const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */ |
1027 | const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */ |
1028 | |
1029 | snew(membed, 1)(membed) = save_calloc("membed", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1029, (1), sizeof(*(membed))); |
1030 | snew(ins_at, 1)(ins_at) = save_calloc("ins_at", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1030, (1), sizeof(*(ins_at))); |
1031 | snew(pos_ins, 1)(pos_ins) = save_calloc("pos_ins", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1031, (1), sizeof(*(pos_ins))); |
1032 | |
1033 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
1034 | { |
1035 | /* get input data out membed file */ |
1036 | membed_input = opt2fn("-membed", nfile, fnm); |
1037 | get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm, |
1038 | &maxwarn, &pieces, &bALLOW_ASYMMETRY); |
1039 | |
1040 | tpr_version = get_tpr_version(ftp2fn(efTPX, nfile, fnm)); |
1041 | if (tpr_version < membed_version) |
1042 | { |
1043 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1043, "Version of *.tpr file to old (%d). " |
1044 | "Rerun grompp with GROMACS version 4.0.3 or newer.\n", tpr_version); |
1045 | } |
1046 | |
1047 | if (!EI_DYNAMICS(inputrec->eI)(((inputrec->eI) == eiMD || ((inputrec->eI) == eiVV || ( inputrec->eI) == eiVVAK)) || ((inputrec->eI) == eiSD1 || (inputrec->eI) == eiSD2) || (inputrec->eI) == eiBD) ) |
1048 | { |
1049 | gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).")_gmx_error("input", "Change integrator to a dynamics integrator in mdp file (e.g. md or sd)." , "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c", 1049); |
1050 | } |
1051 | |
1052 | if (PAR(cr)((cr)->nnodes > 1)) |
1053 | { |
1054 | gmx_input("Sorry, parallel g_membed is not yet fully functional.")_gmx_error("input", "Sorry, parallel g_membed is not yet fully functional." , "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c", 1054); |
1055 | } |
1056 | |
1057 | if (*cpt >= 0) |
1058 | { |
1059 | fprintf(stderrstderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n"); |
1060 | *cpt = -1; |
1061 | } |
1062 | groups = &(mtop->groups); |
1063 | snew(gnames, groups->ngrpname)(gnames) = save_calloc("gnames", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1063, (groups->ngrpname), sizeof(*(gnames))); |
1064 | for (i = 0; i < groups->ngrpname; i++) |
1065 | { |
1066 | gnames[i] = *(groups->grpname[i]); |
1067 | } |
1068 | |
1069 | atoms = gmx_mtop_global_atoms(mtop); |
1070 | snew(mem_p, 1)(mem_p) = save_calloc("mem_p", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1070, (1), sizeof(*(mem_p))); |
1071 | fprintf(stderrstderr, "\nSelect a group to embed in the membrane:\n"); |
1072 | get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins); |
1073 | ins_grp_id = search_string(ins, groups->ngrpname, gnames); |
1074 | fprintf(stderrstderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins); |
1075 | get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name)); |
1076 | |
1077 | pos_ins->pieces = pieces; |
1078 | snew(pos_ins->nidx, pieces)(pos_ins->nidx) = save_calloc("pos_ins->nidx", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1078, (pieces), sizeof(*(pos_ins->nidx))); |
1079 | snew(pos_ins->subindex, pieces)(pos_ins->subindex) = save_calloc("pos_ins->subindex", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1079, (pieces), sizeof(*(pos_ins->subindex))); |
1080 | snew(piecename, pieces)(piecename) = save_calloc("piecename", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1080, (pieces), sizeof(*(piecename))); |
1081 | if (pieces > 1) |
1082 | { |
1083 | fprintf(stderrstderr, "\nSelect pieces to embed:\n"); |
1084 | get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename); |
1085 | } |
1086 | else |
1087 | { |
1088 | /*use whole embedded group*/ |
1089 | snew(pos_ins->nidx, 1)(pos_ins->nidx) = save_calloc("pos_ins->nidx", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1089, (1), sizeof(*(pos_ins->nidx))); |
1090 | snew(pos_ins->subindex, 1)(pos_ins->subindex) = save_calloc("pos_ins->subindex", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1090, (1), sizeof(*(pos_ins->subindex))); |
1091 | pos_ins->nidx[0] = ins_at->nr; |
1092 | pos_ins->subindex[0] = ins_at->index; |
1093 | } |
1094 | |
1095 | if (probe_rad < min_probe_rad) |
1096 | { |
1097 | warn++; |
1098 | fprintf(stderrstderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result " |
1099 | "in overlap between waters and the group to embed, which will result " |
1100 | "in Lincs errors etc.\n\n", warn); |
1101 | } |
1102 | |
1103 | if (xy_fac < min_xy_init) |
1104 | { |
1105 | warn++; |
1106 | fprintf(stderrstderr, "\nWarning %d:\nThe initial size of %s is probably too smal.\n\n", warn, ins); |
1107 | } |
1108 | |
1109 | if (it_xy < min_it_xy) |
1110 | { |
1111 | warn++; |
1112 | fprintf(stderrstderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)" |
1113 | " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy); |
1114 | } |
1115 | |
1116 | if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) ) |
1117 | { |
1118 | warn++; |
1119 | fprintf(stderrstderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)" |
1120 | " is probably too small.\nIncrease -nz or maxwarn.\n\n", warn, ins, it_z); |
1121 | } |
1122 | |
1123 | if (it_xy+it_z > inputrec->nsteps) |
1124 | { |
1125 | warn++; |
1126 | fprintf(stderrstderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the " |
1127 | "number of steps in the tpr.\n\n", warn); |
1128 | } |
1129 | |
1130 | fr_id = -1; |
1131 | if (inputrec->opts.ngfrz == 1) |
1132 | { |
1133 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1133, "You did not specify \"%s\" as a freezegroup.", ins); |
1134 | } |
1135 | |
1136 | for (i = 0; i < inputrec->opts.ngfrz; i++) |
1137 | { |
1138 | tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i]; |
1139 | if (ins_grp_id == tmp_id) |
1140 | { |
1141 | fr_id = tmp_id; |
1142 | fr_i = i; |
1143 | } |
1144 | } |
1145 | |
1146 | if (fr_id == -1) |
1147 | { |
1148 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1148, "\"%s\" not as freezegroup defined in the mdp-file.", ins); |
1149 | } |
1150 | |
1151 | for (i = 0; i < DIM3; i++) |
1152 | { |
1153 | if (inputrec->opts.nFreeze[fr_i][i] != 1) |
1154 | { |
1155 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1155, "freeze dimensions for %s are not Y Y Y\n", ins); |
1156 | } |
1157 | } |
1158 | |
1159 | ng = groups->grps[egcENER].nr; |
1160 | if (ng == 1) |
1161 | { |
1162 | gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group")_gmx_error("input", "No energy groups defined. This is necessary for energy exclusion in the freeze group" , "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c", 1162); |
1163 | } |
1164 | |
1165 | for (i = 0; i < ng; i++) |
1166 | { |
1167 | for (j = 0; j < ng; j++) |
1168 | { |
1169 | if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL(1<<0)) |
1170 | { |
1171 | bExcl = TRUE1; |
1172 | if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || |
1173 | (groups->grps[egcENER].nm_ind[j] != ins_grp_id) ) |
1174 | { |
1175 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1175, "Energy exclusions \"%s\" and \"%s\" do not match the group " |
1176 | "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]], |
1177 | *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins); |
1178 | } |
1179 | } |
1180 | } |
1181 | } |
1182 | |
1183 | if (!bExcl) |
1184 | { |
1185 | gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "_gmx_error("input", "No energy exclusion groups defined. This is necessary for energy exclusion in " "the freeze group", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1186) |
1186 | "the freeze group")_gmx_error("input", "No energy exclusion groups defined. This is necessary for energy exclusion in " "the freeze group", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1186); |
1187 | } |
1188 | |
1189 | /* Obtain the maximum and minimum coordinates of the group to be embedded */ |
1190 | snew(rest_at, 1)(rest_at) = save_calloc("rest_at", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1190, (1), sizeof(*(rest_at))); |
1191 | ins_nat = init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max); |
1192 | /* Check that moleculetypes in insertion group are not part of the rest of the system */ |
1193 | check_types(ins_at, rest_at, mtop); |
1194 | |
1195 | mem_nat = init_mem_at(mem_p, mtop, state->x, state->box, pos_ins); |
Value stored to 'mem_nat' is never read | |
1196 | |
1197 | prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p); |
1198 | if ( (prot_area > prot_vs_box) && ( (state->box[XX0][XX0]*state->box[YY1][YY1]-state->box[XX0][YY1]*state->box[YY1][XX0]) < box_vs_prot) ) |
1199 | { |
1200 | warn++; |
1201 | fprintf(stderrstderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n" |
1202 | "This might cause pressure problems during the growth phase. Just try with\n" |
1203 | "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n" |
1204 | "compressibility in the mdp-file or use no pressure coupling at all.\n\n", warn); |
1205 | } |
1206 | |
1207 | if (warn > maxwarn) |
1208 | { |
1209 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1209, "Too many warnings.\n"); |
1210 | } |
1211 | |
1212 | printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area); |
1213 | printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n" |
1214 | "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area); |
1215 | |
1216 | /* Maximum number of lipids to be removed*/ |
1217 | max_lip_rm = (int)(2*prot_area/mem_p->lip_area); |
1218 | printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm); |
1219 | |
1220 | printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n" |
1221 | "This resizing will be done with respect to the geometrical center of all protein atoms\n" |
1222 | "that span the membrane region, i.e. z between %.3f and %.3f\n\n", |
1223 | xy_fac, z_fac, mem_p->zmin, mem_p->zmax); |
1224 | |
1225 | /* resize the protein by xy and by z if necessary*/ |
1226 | snew(r_ins, ins_at->nr)(r_ins) = save_calloc("r_ins", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1226, (ins_at->nr), sizeof(*(r_ins))); |
1227 | init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY); |
1228 | membed->fac[0] = membed->fac[1] = xy_fac; |
1229 | membed->fac[2] = z_fac; |
1230 | |
1231 | membed->xy_step = (xy_max-xy_fac)/(double)(it_xy); |
1232 | membed->z_step = (z_max-z_fac)/(double)(it_z-1); |
1233 | |
1234 | resize(r_ins, state->x, pos_ins, membed->fac); |
1235 | |
1236 | /* remove overlapping lipids and water from the membrane box*/ |
1237 | /*mark molecules to be removed*/ |
1238 | snew(pbc, 1)(pbc) = save_calloc("pbc", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1238, (1), sizeof(*(pbc))); |
1239 | set_pbc(pbc, inputrec->ePBC, state->box); |
1240 | |
1241 | snew(rm_p, 1)(rm_p) = save_calloc("rm_p", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1241, (1), sizeof(*(rm_p))); |
1242 | lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, mem_p, pos_ins, |
1243 | probe_rad, low_up_rm, bALLOW_ASYMMETRY); |
1244 | lip_rm -= low_up_rm; |
1245 | |
1246 | if (fplog) |
1247 | { |
1248 | for (i = 0; i < rm_p->nr; i++) |
1249 | { |
1250 | fprintf(fplog, "rm mol %d\n", rm_p->mol[i]); |
1251 | } |
1252 | } |
1253 | |
1254 | for (i = 0; i < mtop->nmolblock; i++) |
1255 | { |
1256 | ntype = 0; |
1257 | for (j = 0; j < rm_p->nr; j++) |
1258 | { |
1259 | if (rm_p->block[j] == i) |
1260 | { |
1261 | ntype++; |
1262 | } |
1263 | } |
1264 | printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name)); |
1265 | } |
1266 | |
1267 | if (lip_rm > max_lip_rm) |
1268 | { |
1269 | warn++; |
1270 | fprintf(stderrstderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated " |
1271 | "protein area\nTry making the -xyinit resize factor smaller or increase " |
1272 | "maxwarn.\n\n", warn); |
1273 | } |
1274 | |
1275 | /*remove all lipids and waters overlapping and update all important structures*/ |
1276 | rm_group(groups, mtop, rm_p, state, ins_at, pos_ins); |
1277 | |
1278 | rm_bonded_at = rm_bonded(ins_at, mtop); |
1279 | if (rm_bonded_at != ins_at->nr) |
1280 | { |
1281 | fprintf(stderrstderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, " |
1282 | "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same" |
1283 | "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr); |
1284 | } |
1285 | |
1286 | if (warn > maxwarn) |
1287 | { |
1288 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1288, "Too many warnings.\nIf you are sure these warnings are harmless, " |
1289 | "you can increase -maxwarn"); |
1290 | } |
1291 | |
1292 | if (ftp2bSet(efTOP, nfile, fnm)) |
1293 | { |
1294 | top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop); |
1295 | } |
1296 | |
1297 | sfree(pbc)save_free("pbc", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1297, (pbc)); |
1298 | sfree(rest_at)save_free("rest_at", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1298, (rest_at)); |
1299 | if (pieces > 1) |
1300 | { |
1301 | sfree(piecename)save_free("piecename", "/home/alexxy/Develop/gromacs/src/programs/mdrun/membed.c" , 1301, (piecename)); |
1302 | } |
1303 | |
1304 | membed->it_xy = it_xy; |
1305 | membed->it_z = it_z; |
1306 | membed->pos_ins = pos_ins; |
1307 | membed->r_ins = r_ins; |
1308 | } |
1309 | |
1310 | return membed; |
1311 | } |