4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
27 * For more info, check our website at http://www.gromacs.org
30 * Gromacs Runs One Microsecond At Cannonball Speeds
32 static char *SRCID_genbox_c
= "$Id$";
58 void print_stat(rvec
*x
,int natoms
,matrix box
)
62 for(m
=0;(m
<DIM
);m
++) {
66 for(i
=0;(i
<natoms
);i
++) {
68 xmin
[m
]=min(xmin
[m
],x
[i
][m
]);
69 xmax
[m
]=max(xmax
[m
],x
[i
][m
]);
73 fprintf(stderr
,"DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
74 m
,xmin
[m
],xmax
[m
],box
[m
][m
]);
78 static bool in_box(matrix box
,rvec x
)
80 return( x
[XX
]>=0.0 && x
[XX
]<=box
[XX
][XX
] &&
81 x
[YY
]>=0.0 && x
[YY
]<=box
[YY
][YY
] &&
82 x
[ZZ
]>=0.0 && x
[ZZ
]<=box
[ZZ
][ZZ
] );
85 static void mk_vdw(t_atoms
*a
, real rvdw
[], real r_distance
)
89 /* initialise van der waals arrays of configuration */
90 fprintf(stderr
,"Initialising van der waals distances...\n");
91 for(i
=0; (i
< a
->nr
); i
++)
93 *(a
->resname
[a
->atom
[i
].resnr
]), *(a
->atomname
[i
]),
105 void sort_molecule(t_atoms
**atoms_solvt
,rvec
*x
,rvec
*v
,real
*r
)
107 int atnr
,i
,j
,moltp
=0,nrmoltypes
,resnr
;
108 t_moltypes
*moltypes
;
110 t_atoms
*atoms
,*newatoms
;
111 rvec
*newx
, *newv
=NULL
;
114 fprintf(stderr
,"Sorting configuration\n");
116 atoms
= *atoms_solvt
;
118 /* copy each residue from *atoms to a molecule in *molecule */
123 for (i
=0; i
<atoms
->nr
; i
++) {
124 if ( (i
==0) || (atoms
->atom
[i
].resnr
!= atoms
->atom
[i
-1].resnr
) ) {
125 /* see if this was a molecule type we haven't had yet: */
127 for (j
=0; (j
<nrmoltypes
) && (moltp
==NOTSET
); j
++)
128 if (strcmp(*(atoms
->resname
[atoms
->atom
[i
].resnr
]),
129 moltypes
[j
].name
)==0)
134 srenew(moltypes
,nrmoltypes
);
135 moltypes
[moltp
].name
=*(atoms
->resname
[atoms
->atom
[i
].resnr
]);
137 while ((i
+atnr
<atoms
->nr
) &&
138 (atoms
->atom
[i
].resnr
== atoms
->atom
[i
+atnr
].resnr
))
140 moltypes
[moltp
].natoms
=atnr
;
141 moltypes
[moltp
].nmol
=0;
143 moltypes
[moltp
].nmol
++;
148 fprintf(stderr
,"Found %d%s molecule type%s:\n",
149 nrmoltypes
,nrmoltypes
==1?"":" different",nrmoltypes
==1?"":"s");
150 for(j
=0; j
<nrmoltypes
; j
++) {
152 moltypes
[j
].res0
= 0;
154 moltypes
[j
].res0
= moltypes
[j
-1].res0
+moltypes
[j
-1].nmol
;
155 fprintf(stderr
,"%7s (%4d atoms): %5d residues\n",
156 moltypes
[j
].name
,moltypes
[j
].natoms
,moltypes
[j
].nmol
);
159 /* if we have only 1 moleculetype, we don't have to sort */
161 /* find out which molecules should go where: */
162 moltypes
[0].i
= moltypes
[0].i0
= 0;
163 for(j
=1; j
<nrmoltypes
; j
++) {
166 moltypes
[j
-1].i0
+moltypes
[j
-1].natoms
*moltypes
[j
-1].nmol
;
169 /* now put them there: */
171 init_t_atoms(newatoms
,atoms
->nr
,FALSE
);
172 newatoms
->nres
=atoms
->nres
;
173 snew(newatoms
->resname
,atoms
->nres
);
174 snew(newx
,atoms
->nr
);
175 if (v
) snew(newv
,atoms
->nr
);
176 snew(newr
,atoms
->nr
);
178 for (i
=0; i
<atoms
->nr
; i
++) {
179 resnr
= moltypes
[tps
[i
]].res0
+
180 (moltypes
[tps
[i
]].i
-moltypes
[tps
[i
]].i0
) / moltypes
[tps
[i
]].natoms
;
181 newatoms
->resname
[resnr
] = atoms
->resname
[atoms
->atom
[i
].resnr
];
182 newatoms
->atomname
[moltypes
[tps
[i
]].i
] = atoms
->atomname
[i
];
183 newatoms
->atom
[moltypes
[tps
[i
]].i
] = atoms
->atom
[i
];
184 newatoms
->atom
[moltypes
[tps
[i
]].i
].resnr
= resnr
;
185 copy_rvec(x
[i
],newx
[moltypes
[tps
[i
]].i
]);
186 if (v
) copy_rvec(v
[i
],newv
[moltypes
[tps
[i
]].i
]);
187 newr
[moltypes
[tps
[i
]].i
] = r
[i
];
188 moltypes
[tps
[i
]].i
++;
191 /* put them back into the original arrays and throw away temporary arrays */
192 sfree(atoms
->atomname
);
193 sfree(atoms
->resname
);
196 *atoms_solvt
= newatoms
;
197 for (i
=0; i
<(*atoms_solvt
)->nr
; i
++) {
198 copy_rvec(newx
[i
],x
[i
]);
199 if (v
) copy_rvec(newv
[i
],v
[i
]);
209 void rm_res_pbc(t_atoms
*atoms
, rvec
*x
, matrix box
)
217 for(n
=0; n
<atoms
->nr
; n
++) {
218 if (!is_hydrogen(*(atoms
->atomname
[n
]))) {
222 if ( (n
+1 == atoms
->nr
) ||
223 (atoms
->atom
[n
+1].resnr
!= atoms
->atom
[n
].resnr
) ) {
224 /* if nat==0 we have only hydrogens in the solvent,
225 we take last coordinate as cg */
230 svmul(1.0/nat
,xcg
,xcg
);
231 for(d
=0; d
<DIM
; d
++) {
233 for(i
=start
; i
<=n
; i
++)
237 while (xcg
[d
]>=box
[d
][d
]) {
238 for(i
=start
; i
<=n
; i
++)
250 char *insert_mols(char *mol_insrt
,int nmol_insrt
,int ntry
,int seed
,
251 t_atoms
*atoms
,rvec
**x
,
252 real
**r
,matrix box
, real r_distance
,real rshell
)
254 static char *title_insrt
;
260 real alfa
,beta
,gamma
;
264 /* read number of atoms of insert molecules */
265 get_stx_coordnum(mol_insrt
,&atoms_insrt
.nr
);
266 if (atoms_insrt
.nr
== 0)
267 fatal_error(0,"No molecule in %s, please check your input\n",mol_insrt
);
268 /* allocate memory for atom coordinates of insert molecules */
269 snew(x_insrt
,atoms_insrt
.nr
);
270 snew(r_insrt
,atoms_insrt
.nr
);
271 snew(atoms_insrt
.resname
,atoms_insrt
.nr
);
272 snew(atoms_insrt
.atomname
,atoms_insrt
.nr
);
273 snew(atoms_insrt
.atom
,atoms_insrt
.nr
);
274 atoms_insrt
.pdbinfo
= NULL
;
275 snew(x_n
,atoms_insrt
.nr
);
276 snew(title_insrt
,STRLEN
);
278 /* read residue number, residue names, atomnames, coordinates etc. */
279 fprintf(stderr
,"Reading molecule configuration \n");
280 read_stx_conf(mol_insrt
,title_insrt
,&atoms_insrt
,x_insrt
,NULL
,box_insrt
);
281 fprintf(stderr
,"%s\nContaining %d atoms in %d residue\n",
282 title_insrt
,atoms_insrt
.nr
,atoms_insrt
.nres
);
283 srenew(atoms_insrt
.resname
,atoms_insrt
.nres
);
285 /* initialise van der waals arrays of insert molecules */
286 mk_vdw(&atoms_insrt
, r_insrt
, r_distance
);
288 srenew(atoms
->resname
,(atoms
->nres
+nmol_insrt
));
289 srenew(atoms
->atomname
,(atoms
->nr
+atoms_insrt
.nr
*nmol_insrt
));
290 srenew(atoms
->atom
,(atoms
->nr
+atoms_insrt
.nr
*nmol_insrt
));
291 srenew(*x
,(atoms
->nr
+atoms_insrt
.nr
*nmol_insrt
));
292 srenew(*r
,(atoms
->nr
+atoms_insrt
.nr
*nmol_insrt
));
295 while ((mol
< nmol_insrt
) && (try < ntry
*nmol_insrt
)) {
296 fprintf(stderr
,"\rTry %d",try++);
297 for (i
=0;(i
<atoms_insrt
.nr
);i
++) {
298 if (atoms_insrt
.atom
[i
].resnr
!=0)
299 fatal_error(0,"more then one residue in insert molecules\n"
300 "program terminated\n");
301 copy_rvec(x_insrt
[i
],x_n
[i
]);
303 alfa
=2*M_PI
*rando(&seed
);
304 beta
=2*M_PI
*rando(&seed
);
305 gamma
=2*M_PI
*rando(&seed
);
306 rotate_conf(atoms_insrt
.nr
,x_n
,NULL
,alfa
,beta
,gamma
);
307 offset_x
[XX
]=box
[XX
][XX
]*rando(&seed
);
308 offset_x
[YY
]=box
[YY
][YY
]*rando(&seed
);
309 offset_x
[ZZ
]=box
[ZZ
][ZZ
]*rando(&seed
);
310 gen_box(0,atoms_insrt
.nr
,x_n
,box_insrt
,offset_x
,TRUE
);
311 if (!in_box(box
,x_n
[0]) || !in_box(box
,x_n
[atoms_insrt
.nr
-1]))
315 add_conf(atoms
,x
,NULL
,r
,FALSE
,box
,TRUE
,
316 &atoms_insrt
,x_n
,NULL
,r_insrt
,FALSE
,rshell
);
318 if (atoms
->nr
==(atoms_insrt
.nr
+onr
)) {
320 fprintf(stderr
," success (now %d atoms)!",atoms
->nr
);
323 srenew(atoms
->resname
, atoms
->nres
);
324 srenew(atoms
->atomname
, atoms
->nr
);
325 srenew(atoms
->atom
, atoms
->nr
);
326 srenew(*x
, atoms
->nr
);
327 srenew(*r
, atoms
->nr
);
329 fprintf(stderr
,"\n");
330 /* print number of molecules added */
331 fprintf(stderr
,"Added %d molecules (out of %d requested) of %s\n",
332 mol
,nmol_insrt
,*atoms_insrt
.resname
[0]);
337 void add_solv(char *fn
,t_atoms
*atoms
,rvec
**x
,rvec
**v
,real
**r
,matrix box
,
338 real r_distance
,int *atoms_added
,int *residues_added
,
343 char filename
[STRLEN
];
344 char title_solvt
[STRLEN
];
345 t_atoms
*atoms_solvt
;
346 rvec
*x_solvt
,*v_solvt
=NULL
;
351 strncpy(filename
,libfn(fn
),STRLEN
);
353 get_stx_coordnum(filename
,&(atoms_solvt
->nr
));
354 if (atoms_solvt
->nr
== 0)
355 fatal_error(0,"No solvent in %s, please check your input\n",filename
);
356 snew(x_solvt
,atoms_solvt
->nr
);
357 if (v
) snew(v_solvt
,atoms_solvt
->nr
);
358 snew(r_solvt
,atoms_solvt
->nr
);
359 snew(atoms_solvt
->resname
,atoms_solvt
->nr
);
360 snew(atoms_solvt
->atomname
,atoms_solvt
->nr
);
361 snew(atoms_solvt
->atom
,atoms_solvt
->nr
);
362 atoms_solvt
->pdbinfo
= NULL
;
363 fprintf(stderr
,"Reading solvent configuration%s\n",
364 v_solvt
?" and velocities":"");
365 read_stx_conf(filename
,title_solvt
,atoms_solvt
,x_solvt
,v_solvt
,box_solvt
);
366 fprintf(stderr
,"\"%s\"\n",title_solvt
);
367 fprintf(stderr
,"solvent configuration contains %d atoms in %d residues\n",
368 atoms_solvt
->nr
,atoms_solvt
->nres
);
369 fprintf(stderr
,"\n");
371 /* apply pbc for solvent configuration for whole molecules */
372 rm_res_pbc(atoms_solvt
,x_solvt
,box_solvt
);
374 /* initialise van der waals arrays of solvent configuration */
375 mk_vdw(atoms_solvt
,r_solvt
,r_distance
);
377 /* calculate the box multiplication factors n_box[0...DIM] */
379 for (i
=0; (i
< DIM
);i
++) {
381 while (n_box
[i
]*box_solvt
[i
][i
] < box
[i
][i
])
385 fprintf(stderr
,"Will generate new solvent configuration of %dx%dx%d boxes\n",
386 n_box
[XX
],n_box
[YY
],n_box
[ZZ
]);
388 /* realloc atoms_solvt for the new solvent configuration */
389 srenew(atoms_solvt
->resname
,atoms_solvt
->nres
*nmol
);
390 srenew(atoms_solvt
->atomname
,atoms_solvt
->nr
*nmol
);
391 srenew(atoms_solvt
->atom
,atoms_solvt
->nr
*nmol
);
392 srenew(x_solvt
,atoms_solvt
->nr
*nmol
);
393 if (v_solvt
) srenew(v_solvt
,atoms_solvt
->nr
*nmol
);
394 srenew(r_solvt
,atoms_solvt
->nr
*nmol
);
396 /* generate a new solvent configuration */
397 genconf(atoms_solvt
,x_solvt
,v_solvt
,r_solvt
,box_solvt
,n_box
);
400 print_stat(x_solvt
,atoms_solvt
->nr
,box_solvt
);
404 print_stat(x_solvt
,atoms_solvt
->nr
,box_solvt
);
406 /* Sort the solvent mixture, not the protein... */
407 sort_molecule(&atoms_solvt
,x_solvt
,v_solvt
,r_solvt
);
409 /* add the two configurations */
412 add_conf(atoms
,x
,v
,r
,TRUE
,box
,FALSE
,
413 atoms_solvt
,x_solvt
,v_solvt
,r_solvt
,TRUE
,rshell
);
414 *atoms_added
=atoms
->nr
-onr
;
415 *residues_added
=atoms
->nres
-onres
;
420 fprintf(stderr
,"Generated solvent containing %d atoms in %d residues\n",
421 *atoms_added
,*residues_added
);
424 char *read_prot(char *confin
,t_atoms
*atoms
,rvec
**x
,rvec
**v
,real
**r
,
425 matrix box
,real r_distance
)
431 get_stx_coordnum(confin
,&natoms
);
433 /* allocate memory for atom coordinates of configuration 1 */
435 if (v
) snew(*v
,natoms
);
437 init_t_atoms(atoms
,natoms
,FALSE
);
439 /* read residue number, residue names, atomnames, coordinates etc. */
440 fprintf(stderr
,"Reading solute configuration%s\n",v
?" and velocities":"");
441 read_stx_conf(confin
,title
,atoms
,*x
,v
?*v
:NULL
,box
);
442 fprintf(stderr
,"%s\nContaining %d atoms in %d residues\n",
443 title
,atoms
->nr
,atoms
->nres
);
445 /* initialise van der waals arrays of configuration 1 */
446 mk_vdw(atoms
,*r
,r_distance
);
451 static void update_top(t_atoms
*atoms
,matrix box
,int NFILE
,t_filenm fnm
[])
453 #define TEMP_FILENM "temp.top"
455 char buf
[STRLEN
],buf2
[STRLEN
],*temp
,*topinout
;
457 bool bSystem
,bMolecules
,bSkip
;
462 for(i
=0; (i
<atoms
->nres
); i
++) {
463 /* calculate number of SOLvent molecules */
464 if ( (strcmp(*atoms
->resname
[i
],"SOL")==0) ||
465 (strcmp(*atoms
->resname
[i
],"WAT")==0) ||
466 (strcmp(*atoms
->resname
[i
],"HOH")==0) )
470 for(i
=0; (i
<atoms
->nr
); i
++)
471 mtot
+= get_mass(*atoms
->resname
[atoms
->atom
[i
].resnr
],
472 *atoms
->atomname
[i
]);
476 fprintf(stderr
,"Volume : %10g (nm^3)\n",vol
);
477 fprintf(stderr
,"Density : %10g (g/l)\n",
478 (mtot
*1e24
)/(AVOGADRO
*vol
));
479 fprintf(stderr
,"Number of SOL molecules: %5d \n\n",nsol
);
481 /* open topology file and append sol molecules */
482 topinout
= ftp2fn(efTOP
,NFILE
,fnm
);
483 if ( ftp2bSet(efTOP
,NFILE
,fnm
) ) {
484 fprintf(stderr
,"Processing topology\n");
485 fpin
= ffopen(topinout
,"r");
486 fpout
= ffopen(TEMP_FILENM
,"w");
488 bSystem
= bMolecules
= FALSE
;
489 while (fgets(buf
, STRLEN
, fpin
)) {
493 if ((temp
=strchr(buf2
,'\n')) != NULL
)
498 if ((temp
=strchr(buf2
,'\n')) != NULL
)
501 if (buf2
[strlen(buf2
)-1]==']') {
502 buf2
[strlen(buf2
)-1]='\0';
505 bSystem
=(strcasecmp(buf2
,"system")==0);
506 bMolecules
=(strcasecmp(buf2
,"molecules")==0);
508 } else if (bSystem
&& nsol
&& (buf
[0]!=';') ) {
509 /* if sol present, append "in water" to system name */
511 if (buf2
[0] && (!strstr(buf2
," water")) ) {
512 sprintf(buf
,"%s in water\n",buf2
);
515 } else if (bMolecules
) {
516 /* check if this is a line with solvent molecules */
517 sscanf(buf
,"%s",buf2
);
518 if (strcmp(buf2
,"SOL")==0) {
519 sscanf(buf
,"%*s %d",&i
);
528 if ((temp
=strchr(buf
,'\n')) != NULL
)
530 fprintf(stdout
,"Removing line #%d '%s' from topology file (%s)\n",
533 fprintf(fpout
,"%s",buf
);
537 fprintf(stdout
,"Adding line for %d solute molecules to "
538 "topology file (%s)\n",nsol
,topinout
);
539 fprintf(fpout
,"%-15s %5d\n","SOL",nsol
);
542 /* use ffopen to generate backup of topinout */
543 fpout
=ffopen(topinout
,"w");
545 rename(TEMP_FILENM
,topinout
);
550 int main(int argc
,char *argv
[])
552 static char *desc
[] = {
553 "Genbox can do one of 3 things:[PAR]",
555 "1) Generate a box of solvent. Specify -cs and -box. Or specify -cs and",
556 "-cp with a structure file with a box, but without atoms.[PAR]",
558 "2) Solvate a solute configuration, eg. a protein, in a bath of solvent ",
559 "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
560 "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
561 "unless [TT]-box[tt] is set, which also centers the solute.",
562 "The program [TT]editconf[tt] has more sophisticated options to change",
563 "the box and center the solute.",
564 "Solvent molecules are removed from the box where the ",
565 "distance between any atom of the solute molecule(s) and any atom of ",
566 "the solvent molecule is less than the sum of the VanderWaals radii of ",
567 "both atoms. A database ([TT]vdwradii.dat[tt]) of VanderWaals radii is ",
568 "read by the program, atoms not in the database are ",
569 "assigned a default distance [TT]-vdw[tt].[PAR]",
571 "3) Insert a number ([TT]-nmol[tt]) of extra molecules ([TT]-ci[tt]) ",
572 "at random positions.",
573 "The program iterates until [TT]nmol[tt] molecules",
574 "have been inserted in the box. To test whether an insertion is ",
575 "successful the same VanderWaals criterium is used as for removal of ",
576 "solvent molecules. When no appropriately ",
577 "sized holes (holes that can hold an extra molecule) are available the ",
578 "program tries for [TT]-nmol[tt] * [TT]-try[tt] times before giving up. ",
579 "Increase -try if you have several small holes to fill.[PAR]",
581 "The default solvent is Simple Point Charge water (SPC), with coordinates ",
582 "from [TT]$GMXLIB/spc216.gro[tt]. Other",
583 "solvents are also supported, as well as mixed solvents. The",
584 "only restriction to solvent types is that a solvent molecule consists",
585 "of exactly one residue. The residue information in the coordinate",
586 "files is used, and should therefore be more or less consistent.",
587 "In practice this means that two subsequent solvent molecules in the ",
588 "solvent coordinate file should have different residue number.",
589 "The box of solute is built by stacking the coordinates read from",
590 "the coordinate file. This means that these coordinates should be ",
591 "equlibrated in periodic boundary conditions to ensure a good",
592 "alignment of molecules on the stacking interfaces.[PAR]",
594 "The program can optionally rotate the solute molecule to align the",
595 "longest molecule axis along a box edge. This way the amount of solvent",
596 "molecules necessary is reduced.",
597 "It should be kept in mind that this only works for",
598 "short simulations, as eg. an alpha-helical peptide in solution can ",
599 "rotate over 90 degrees, within 500 ps. In general it is therefore ",
600 "better to make a more or less cubic box.[PAR]",
602 "Setting -shell larger than zero will place a layer of water of",
603 "the specified thickness (nm) around the solute. Hint: it is a good",
604 "idea to put the protein in the center of a box first (using editconf).",
607 "Finally, genbox will optionally remove lines from your topology file in ",
608 "which a number of solvent molecules is already added, and adds a ",
609 "line with the total number of solvent molecules in your coordinate file."
612 static char *bugs
[] = {
613 "Molecules must be whole in the initial configurations.",
614 "At the moment -ci only works when inserting one molecule."
618 bool bSol
,bProt
,bBox
;
619 char *conf_prot
,*confout
;
624 /* protein configuration data */
630 /* other data types */
631 int atoms_added
,residues_added
;
634 { efSTX
, "-cp", "protein", ffOPTRD
},
635 { efSTX
, "-cs", "spc216", ffLIBOPTRD
},
636 { efSTX
, "-ci", "insert", ffOPTRD
},
637 { efSTO
, NULL
, NULL
, ffWRITE
},
638 { efTOP
, NULL
, NULL
, ffOPTRW
},
640 #define NFILE asize(fnm)
642 static int nmol_ins
=0,nmol_try
=10,seed
=1997;
643 static real r_distance
=0.105,r_shell
=0;
644 static rvec new_box
={0.0,0.0,0.0};
645 static bool bReadV
=FALSE
;
647 { "-box", FALSE
, etRVEC
, {new_box
},
649 { "-nmol", FALSE
, etINT
, {&nmol_ins
},
650 "no of extra molecules to insert" },
651 { "-try", FALSE
, etINT
, {&nmol_try
},
652 "try inserting -nmol*-try times" },
653 { "-seed", FALSE
, etINT
, {&seed
},
654 "random generator seed"},
655 { "-vdwd", FALSE
, etREAL
, {&r_distance
},
656 "default vdwaals distance"},
657 { "-shell", FALSE
, etREAL
, {&r_shell
},
658 "thickness of optional water layer around solute" },
659 { "-vel", FALSE
, etBOOL
, {&bReadV
},
660 "HIDDENkeep velocities from input solute and solvent" }
663 CopyRight(stderr
,argv
[0]);
664 parse_common_args(&argc
,argv
, PCA_BE_NICE
,NFILE
,fnm
,asize(pa
),pa
,
665 asize(desc
),desc
,asize(bugs
),bugs
);
667 bInsert
= opt2bSet("-ci",NFILE
,fnm
) && (nmol_ins
> 0);
668 bSol
= opt2bSet("-cs",NFILE
,fnm
);
669 bProt
= opt2bSet("-cp",NFILE
,fnm
);
670 bBox
= opt2parg_bSet("-box",asize(pa
),pa
);
673 if (bInsert
&& nmol_ins
<=0)
674 fatal_error(0,"When specifying inserted molecules (-ci), "
675 "-nmol must be larger than 0");
677 fatal_error(0,"When no solute (-cp) is specified, "
678 "a box size (-box) must be specified");
681 /*generate a solute configuration */
682 conf_prot
= opt2fn("-cp",NFILE
,fnm
);
683 title
= read_prot(conf_prot
,&atoms
,&x
,bReadV
?&v
:NULL
,&r
,box
,r_distance
);
685 fprintf(stderr
,"Note: no velocities found\n");
687 fprintf(stderr
,"Note: no atoms in %s\n",conf_prot
);
703 box
[XX
][XX
]=new_box
[XX
];
704 box
[YY
][YY
]=new_box
[YY
];
705 box
[ZZ
][ZZ
]=new_box
[ZZ
];
708 fatal_error(0,"Undefined solute box.\nCreate one with editconf "
709 "or give explicit -box command line option");
713 /* add nmol_ins molecules of atoms_ins
714 in random orientation at random place */
716 title_ins
= insert_mols(opt2fn("-ci",NFILE
,fnm
),nmol_ins
,nmol_try
,seed
,
717 &atoms
,&x
,&r
,box
,r_distance
,r_shell
);
719 title_ins
= strdup("Generated by genbox");
723 add_solv(opt2fn("-cs",NFILE
,fnm
),&atoms
,&x
,v
?&v
:NULL
,&r
,box
,
724 r_distance
,&atoms_added
,&residues_added
,r_shell
);
726 /* write new configuration 1 to file confout */
727 confout
= ftp2fn(efSTO
,NFILE
,fnm
);
728 fprintf(stderr
,"Writing generated configuration to %s\n",confout
);
730 write_sto_conf(confout
,title
,&atoms
,x
,v
,box
);
731 /* print box sizes and box type to stderr */
732 fprintf(stderr
,"%s\n",title
);
734 write_sto_conf(confout
,title_ins
,&atoms
,x
,v
,box
);
736 /* print size of generated configuration */
737 fprintf(stderr
,"\nOutput configuration contains %d atoms in %d residues\n",
738 atoms
.nr
,atoms
.nres
);
739 update_top(&atoms
,box
,NFILE
,fnm
);