Bug Summary

File:programs/mdrun/membed.c
Location:line 169, column 9
Description:Value stored to 'mol_id' is never read

Annotated Source Code

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 */
60typedef 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 */
71struct 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 */
82typedef 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 */
95typedef 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. */
103static 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 */
124static 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
143static 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. */
158static 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);
Value stored to 'mol_id' is never read
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 */
189static 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
224static 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 */
251static 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 */
341static 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
378static 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
455static 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 */
515static 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. */
534static 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)*/
691static 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 */
821int 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 */
889static 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
971void 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
986gmx_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);
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}