/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
/*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
- *
+ *
* If you want to redistribute modifications, please consider that
* scientific software is very special. Version control is crucial -
* bugs must be traceable. We will be happy to consider code for
* inclusion in the official distribution, but derived work must not
* be called official GROMACS. Details are found in the README & COPYING
* files - if they are missing, get the official version at www.gromacs.org.
- *
+ *
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
- *
+ *
* For more info, check our website at http://www.gromacs.org
- *
+ *
* And Hey:
* Green Red Orange Magenta Azure Cyan Skyblue
*/
#include "xvgr.h"
#ifdef DEBUG
-void print_stat(rvec *x,int natoms,matrix box)
+void print_stat(rvec *x, int natoms, matrix box)
{
- int i,m;
- rvec xmin,xmax;
- for(m=0;(m<DIM);m++) {
- xmin[m]=x[0][m];
- xmax[m]=x[0][m];
- }
- for(i=0;(i<natoms);i++) {
- for (m=0;m<DIM;m++) {
- xmin[m]=min(xmin[m],x[i][m]);
- xmax[m]=max(xmax[m],x[i][m]);
- }
- }
- for(m=0;(m<DIM);m++)
- fprintf(stderr,"DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
- m,xmin[m],xmax[m],box[m][m]);
+ int i, m;
+ rvec xmin, xmax;
+ for (m = 0; (m < DIM); m++)
+ {
+ xmin[m] = x[0][m];
+ xmax[m] = x[0][m];
+ }
+ for (i = 0; (i < natoms); i++)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ xmin[m] = min(xmin[m], x[i][m]);
+ xmax[m] = max(xmax[m], x[i][m]);
+ }
+ }
+ for (m = 0; (m < DIM); m++)
+ {
+ fprintf(stderr, "DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
+ m, xmin[m], xmax[m], box[m][m]);
+ }
}
#endif
-static gmx_bool in_box(t_pbc *pbc,rvec x)
+static gmx_bool in_box(t_pbc *pbc, rvec x)
{
- rvec box_center,dx;
- int shift;
-
- /* pbc_dx_aiuc only works correctly with the rectangular box center */
- calc_box_center(ecenterRECT,pbc->box,box_center);
-
- shift = pbc_dx_aiuc(pbc,x,box_center,dx);
-
- return (shift == CENTRAL);
+ rvec box_center, dx;
+ int shift;
+
+ /* pbc_dx_aiuc only works correctly with the rectangular box center */
+ calc_box_center(ecenterRECT, pbc->box, box_center);
+
+ shift = pbc_dx_aiuc(pbc, x, box_center, dx);
+
+ return (shift == CENTRAL);
}
-static void mk_vdw(t_atoms *a,real rvdw[],gmx_atomprop_t aps,
- real r_distance)
+static void mk_vdw(t_atoms *a, real rvdw[], gmx_atomprop_t aps,
+ real r_distance)
{
- int i;
-
- /* initialise van der waals arrays of configuration */
- fprintf(stderr,"Initialising van der waals distances...\n");
- for(i=0; (i < a->nr); i++)
- if (!gmx_atomprop_query(aps,epropVDW,
- *(a->resinfo[a->atom[i].resind].name),
- *(a->atomname[i]),&(rvdw[i])))
- rvdw[i] = r_distance;
+ int i;
+
+ /* initialise van der waals arrays of configuration */
+ fprintf(stderr, "Initialising van der waals distances...\n");
+ for (i = 0; (i < a->nr); i++)
+ {
+ if (!gmx_atomprop_query(aps, epropVDW,
+ *(a->resinfo[a->atom[i].resind].name),
+ *(a->atomname[i]), &(rvdw[i])))
+ {
+ rvdw[i] = r_distance;
+ }
+ }
}
typedef struct {
- char *name;
- int natoms;
- int nmol;
- int i,i0;
- int res0;
+ char *name;
+ int natoms;
+ int nmol;
+ int i, i0;
+ int res0;
} t_moltypes;
-void sort_molecule(t_atoms **atoms_solvt,rvec *x,rvec *v,real *r)
+void sort_molecule(t_atoms **atoms_solvt, rvec *x, rvec *v, real *r)
{
- int atnr,i,j,moltp=0,nrmoltypes,resi_o,resi_n,resnr;
- t_moltypes *moltypes;
- t_atoms *atoms,*newatoms;
- rvec *newx, *newv=NULL;
- real *newr;
-
- fprintf(stderr,"Sorting configuration\n");
-
- atoms = *atoms_solvt;
-
- /* copy each residue from *atoms to a molecule in *molecule */
- moltypes=NULL;
- nrmoltypes=0;
- for (i=0; i<atoms->nr; i++) {
- if ( (i==0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) ) {
- /* see if this was a molecule type we haven't had yet: */
- moltp=NOTSET;
- for (j=0; (j<nrmoltypes) && (moltp==NOTSET); j++)
- if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
- moltypes[j].name)==0)
- moltp=j;
- if (moltp==NOTSET) {
- moltp=nrmoltypes;
- nrmoltypes++;
- srenew(moltypes,nrmoltypes);
- moltypes[moltp].name=*(atoms->resinfo[atoms->atom[i].resind].name);
- atnr = 0;
- while ((i+atnr<atoms->nr) &&
- (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
- atnr++;
- moltypes[moltp].natoms=atnr;
- moltypes[moltp].nmol=0;
- }
- moltypes[moltp].nmol++;
- }
- }
-
- fprintf(stderr,"Found %d%s molecule type%s:\n",
- nrmoltypes,nrmoltypes==1?"":" different",nrmoltypes==1?"":"s");
- for(j=0; j<nrmoltypes; j++) {
- if (j==0)
- moltypes[j].res0 = 0;
- else
- moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
- fprintf(stderr,"%7s (%4d atoms): %5d residues\n",
- moltypes[j].name,moltypes[j].natoms,moltypes[j].nmol);
- }
-
- /* if we have only 1 moleculetype, we don't have to sort */
- if (nrmoltypes>1) {
- /* find out which molecules should go where: */
- moltypes[0].i = moltypes[0].i0 = 0;
- for(j=1; j<nrmoltypes; j++) {
- moltypes[j].i =
- moltypes[j].i0 =
- moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
+ int atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
+ t_moltypes *moltypes;
+ t_atoms *atoms, *newatoms;
+ rvec *newx, *newv = NULL;
+ real *newr;
+
+ fprintf(stderr, "Sorting configuration\n");
+
+ atoms = *atoms_solvt;
+
+ /* copy each residue from *atoms to a molecule in *molecule */
+ moltypes = NULL;
+ nrmoltypes = 0;
+ for (i = 0; i < atoms->nr; i++)
+ {
+ if ( (i == 0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) )
+ {
+ /* see if this was a molecule type we haven't had yet: */
+ moltp = NOTSET;
+ for (j = 0; (j < nrmoltypes) && (moltp == NOTSET); j++)
+ {
+ if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
+ moltypes[j].name) == 0)
+ {
+ moltp = j;
+ }
+ }
+ if (moltp == NOTSET)
+ {
+ moltp = nrmoltypes;
+ nrmoltypes++;
+ srenew(moltypes, nrmoltypes);
+ moltypes[moltp].name = *(atoms->resinfo[atoms->atom[i].resind].name);
+ atnr = 0;
+ while ((i+atnr < atoms->nr) &&
+ (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
+ {
+ atnr++;
+ }
+ moltypes[moltp].natoms = atnr;
+ moltypes[moltp].nmol = 0;
+ }
+ moltypes[moltp].nmol++;
+ }
}
-
- /* now put them there: */
- snew(newatoms,1);
- init_t_atoms(newatoms,atoms->nr,FALSE);
- newatoms->nres=atoms->nres;
- snew(newatoms->resinfo,atoms->nres);
- snew(newx,atoms->nr);
- if (v) snew(newv,atoms->nr);
- snew(newr,atoms->nr);
-
- resi_n = 0;
- resnr = 1;
- j = 0;
- for(moltp=0; moltp<nrmoltypes; moltp++) {
- i = 0;
- while (i < atoms->nr) {
- resi_o = atoms->atom[i].resind;
- if (strcmp(*atoms->resinfo[resi_o].name,moltypes[moltp].name) == 0) {
- /* Copy the residue info */
- newatoms->resinfo[resi_n] = atoms->resinfo[resi_o];
- newatoms->resinfo[resi_n].nr = resnr;
- /* Copy the atom info */
- do {
- newatoms->atom[j] = atoms->atom[i];
- newatoms->atomname[j] = atoms->atomname[i];
- newatoms->atom[j].resind = resi_n;
- copy_rvec(x[i],newx[j]);
- if (v != NULL) {
- copy_rvec(v[i],newv[j]);
- }
- newr[j] = r[i];
- i++;
- j++;
- } while (i < atoms->nr && atoms->atom[i].resind == resi_o);
- /* Increase the new residue counters */
- resi_n++;
- resnr++;
- } else {
- /* Skip this residue */
- do {
- i++;
- } while (i < atoms->nr && atoms->atom[i].resind == resi_o);
- }
- }
+
+ fprintf(stderr, "Found %d%s molecule type%s:\n",
+ nrmoltypes, nrmoltypes == 1 ? "" : " different", nrmoltypes == 1 ? "" : "s");
+ for (j = 0; j < nrmoltypes; j++)
+ {
+ if (j == 0)
+ {
+ moltypes[j].res0 = 0;
+ }
+ else
+ {
+ moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
+ }
+ fprintf(stderr, "%7s (%4d atoms): %5d residues\n",
+ moltypes[j].name, moltypes[j].natoms, moltypes[j].nmol);
}
-
- /* put them back into the original arrays and throw away temporary arrays */
- sfree(atoms->atomname);
- sfree(atoms->resinfo);
- sfree(atoms->atom);
- sfree(atoms);
- *atoms_solvt = newatoms;
- for (i=0; i<(*atoms_solvt)->nr; i++) {
- copy_rvec(newx[i],x[i]);
- if (v) copy_rvec(newv[i],v[i]);
- r[i]=newr[i];
+
+ /* if we have only 1 moleculetype, we don't have to sort */
+ if (nrmoltypes > 1)
+ {
+ /* find out which molecules should go where: */
+ moltypes[0].i = moltypes[0].i0 = 0;
+ for (j = 1; j < nrmoltypes; j++)
+ {
+ moltypes[j].i =
+ moltypes[j].i0 =
+ moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
+ }
+
+ /* now put them there: */
+ snew(newatoms, 1);
+ init_t_atoms(newatoms, atoms->nr, FALSE);
+ newatoms->nres = atoms->nres;
+ snew(newatoms->resinfo, atoms->nres);
+ snew(newx, atoms->nr);
+ if (v)
+ {
+ snew(newv, atoms->nr);
+ }
+ snew(newr, atoms->nr);
+
+ resi_n = 0;
+ resnr = 1;
+ j = 0;
+ for (moltp = 0; moltp < nrmoltypes; moltp++)
+ {
+ i = 0;
+ while (i < atoms->nr)
+ {
+ resi_o = atoms->atom[i].resind;
+ if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
+ {
+ /* Copy the residue info */
+ newatoms->resinfo[resi_n] = atoms->resinfo[resi_o];
+ newatoms->resinfo[resi_n].nr = resnr;
+ /* Copy the atom info */
+ do
+ {
+ newatoms->atom[j] = atoms->atom[i];
+ newatoms->atomname[j] = atoms->atomname[i];
+ newatoms->atom[j].resind = resi_n;
+ copy_rvec(x[i], newx[j]);
+ if (v != NULL)
+ {
+ copy_rvec(v[i], newv[j]);
+ }
+ newr[j] = r[i];
+ i++;
+ j++;
+ }
+ while (i < atoms->nr && atoms->atom[i].resind == resi_o);
+ /* Increase the new residue counters */
+ resi_n++;
+ resnr++;
+ }
+ else
+ {
+ /* Skip this residue */
+ do
+ {
+ i++;
+ }
+ while (i < atoms->nr && atoms->atom[i].resind == resi_o);
+ }
+ }
+ }
+
+ /* put them back into the original arrays and throw away temporary arrays */
+ sfree(atoms->atomname);
+ sfree(atoms->resinfo);
+ sfree(atoms->atom);
+ sfree(atoms);
+ *atoms_solvt = newatoms;
+ for (i = 0; i < (*atoms_solvt)->nr; i++)
+ {
+ copy_rvec(newx[i], x[i]);
+ if (v)
+ {
+ copy_rvec(newv[i], v[i]);
+ }
+ r[i] = newr[i];
+ }
+ sfree(newx);
+ if (v)
+ {
+ sfree(newv);
+ }
+ sfree(newr);
}
- sfree(newx);
- if (v) sfree(newv);
- sfree(newr);
- }
- sfree(moltypes);
+ sfree(moltypes);
}
void rm_res_pbc(t_atoms *atoms, rvec *x, matrix box)
{
- int i,start,n,d,nat;
- rvec xcg;
-
- start=0;
- nat=0;
- clear_rvec(xcg);
- for(n=0; n<atoms->nr; n++) {
- if (!is_hydrogen(*(atoms->atomname[n]))) {
- nat++;
- rvec_inc(xcg,x[n]);
- }
- if ( (n+1 == atoms->nr) ||
- (atoms->atom[n+1].resind != atoms->atom[n].resind) ) {
- /* if nat==0 we have only hydrogens in the solvent,
- we take last coordinate as cg */
- if (nat==0) {
- nat=1;
- copy_rvec(x[n],xcg);
- }
- svmul(1.0/nat,xcg,xcg);
- for(d=0; d<DIM; d++) {
- while (xcg[d]<0) {
- for(i=start; i<=n; i++)
- x[i][d]+=box[d][d];
- xcg[d]+=box[d][d];
- }
- while (xcg[d]>=box[d][d]) {
- for(i=start; i<=n; i++)
- x[i][d]-=box[d][d];
- xcg[d]-=box[d][d];
- }
- }
- start=n+1;
- nat=0;
- clear_rvec(xcg);
+ int i, start, n, d, nat;
+ rvec xcg;
+
+ start = 0;
+ nat = 0;
+ clear_rvec(xcg);
+ for (n = 0; n < atoms->nr; n++)
+ {
+ if (!is_hydrogen(*(atoms->atomname[n])))
+ {
+ nat++;
+ rvec_inc(xcg, x[n]);
+ }
+ if ( (n+1 == atoms->nr) ||
+ (atoms->atom[n+1].resind != atoms->atom[n].resind) )
+ {
+ /* if nat==0 we have only hydrogens in the solvent,
+ we take last coordinate as cg */
+ if (nat == 0)
+ {
+ nat = 1;
+ copy_rvec(x[n], xcg);
+ }
+ svmul(1.0/nat, xcg, xcg);
+ for (d = 0; d < DIM; d++)
+ {
+ while (xcg[d] < 0)
+ {
+ for (i = start; i <= n; i++)
+ {
+ x[i][d] += box[d][d];
+ }
+ xcg[d] += box[d][d];
+ }
+ while (xcg[d] >= box[d][d])
+ {
+ for (i = start; i <= n; i++)
+ {
+ x[i][d] -= box[d][d];
+ }
+ xcg[d] -= box[d][d];
+ }
+ }
+ start = n+1;
+ nat = 0;
+ clear_rvec(xcg);
+ }
}
- }
}
/* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
static gmx_bool
allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
int ePBC, matrix box,
- t_atoms *atoms_insrt, rvec *x_n,real *r_insrt)
+ t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
{
- int i,j;
- rvec dx;
- real n2,r2;
+ int i, j;
+ rvec dx;
+ real n2, r2;
t_pbc pbc;
- set_pbc(&pbc,ePBC,box);
- for (i=0; i<atoms->nr; i++)
+ set_pbc(&pbc, ePBC, box);
+ for (i = 0; i < atoms->nr; i++)
{
- for (j=0; j<atoms_insrt->nr; j++)
+ for (j = 0; j < atoms_insrt->nr; j++)
{
- pbc_dx(&pbc,x[i],x_n[j],dx);
+ pbc_dx(&pbc, x[i], x_n[j], dx);
n2 = norm2(dx);
r2 = sqr(r[i]+r_insrt[j]);
- if (n2<r2)
+ if (n2 < r2)
+ {
return FALSE;
+ }
}
}
return TRUE;
}
/* enum for random rotations of inserted solutes */
-enum {en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR};
+enum {
+ en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR
+};
-static char *insert_mols(const char *mol_insrt,int nmol_insrt,int ntry,int seed,
- t_atoms *atoms,rvec **x,real **r,int ePBC,matrix box,
- gmx_atomprop_t aps,real r_distance,real rshell,
+static char *insert_mols(const char *mol_insrt, int nmol_insrt, int ntry, int seed,
+ t_atoms *atoms, rvec **x, real **r, int ePBC, matrix box,
+ gmx_atomprop_t aps, real r_distance, real rshell,
const output_env_t oenv,
- const char* posfn, const rvec deltaR,int enum_rot,
+ const char* posfn, const rvec deltaR, int enum_rot,
gmx_bool bCheckAllPairDist)
{
- t_pbc pbc;
- static char *title_insrt;
- t_atoms atoms_insrt;
- rvec *x_insrt,*x_n;
- real *r_insrt;
- int ePBC_insrt;
- matrix box_insrt;
- int i,mol,onr,ncol;
- real alfa=0.,beta=0.,gamma=0.;
- rvec offset_x;
- int trial;
- double **rpos;
-
- set_pbc(&pbc,ePBC,box);
-
- /* read number of atoms of insert molecules */
- get_stx_coordnum(mol_insrt,&atoms_insrt.nr);
- if (atoms_insrt.nr == 0)
- gmx_fatal(FARGS,"No molecule in %s, please check your input\n",mol_insrt);
- /* allocate memory for atom coordinates of insert molecules */
- snew(x_insrt,atoms_insrt.nr);
- snew(r_insrt,atoms_insrt.nr);
- snew(atoms_insrt.resinfo,atoms_insrt.nr);
- snew(atoms_insrt.atomname,atoms_insrt.nr);
- snew(atoms_insrt.atom,atoms_insrt.nr);
- atoms_insrt.pdbinfo = NULL;
- snew(x_n,atoms_insrt.nr);
- snew(title_insrt,STRLEN);
-
- /* read residue number, residue names, atomnames, coordinates etc. */
- fprintf(stderr,"Reading molecule configuration \n");
- read_stx_conf(mol_insrt,title_insrt,&atoms_insrt,x_insrt,NULL,
- &ePBC_insrt,box_insrt);
- fprintf(stderr,"%s\nContaining %d atoms in %d residue\n",
- title_insrt,atoms_insrt.nr,atoms_insrt.nres);
- srenew(atoms_insrt.resinfo,atoms_insrt.nres);
-
- /* initialise van der waals arrays of insert molecules */
- mk_vdw(&atoms_insrt,r_insrt,aps,r_distance);
-
- /* With -ip, take nmol_insrt from file posfn */
- if (posfn != NULL)
- {
- nmol_insrt=read_xvg(posfn,&rpos,&ncol);
- if (ncol != 3)
- gmx_fatal(FARGS,"Expected 3 columns (x/y/z coordinates) in file %s\n",ncol,posfn);
- fprintf(stderr,"Read %d positions from file %s\n\n",nmol_insrt,posfn);
- }
-
- srenew(atoms->resinfo,(atoms->nres+nmol_insrt*atoms_insrt.nres));
- srenew(atoms->atomname,(atoms->nr+atoms_insrt.nr*nmol_insrt));
- srenew(atoms->atom,(atoms->nr+atoms_insrt.nr*nmol_insrt));
- srenew(*x,(atoms->nr+atoms_insrt.nr*nmol_insrt));
- srenew(*r,(atoms->nr+atoms_insrt.nr*nmol_insrt));
-
- trial=mol=0;
- while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt)) {
- fprintf(stderr,"\rTry %d",trial++);
- for (i=0;(i<atoms_insrt.nr);i++) {
- copy_rvec(x_insrt[i],x_n[i]);
- }
- switch (enum_rot)
+ t_pbc pbc;
+ static char *title_insrt;
+ t_atoms atoms_insrt;
+ rvec *x_insrt, *x_n;
+ real *r_insrt;
+ int ePBC_insrt;
+ matrix box_insrt;
+ int i, mol, onr, ncol;
+ real alfa = 0., beta = 0., gamma = 0.;
+ rvec offset_x;
+ int trial;
+ double **rpos;
+
+ set_pbc(&pbc, ePBC, box);
+
+ /* read number of atoms of insert molecules */
+ get_stx_coordnum(mol_insrt, &atoms_insrt.nr);
+ if (atoms_insrt.nr == 0)
{
- case en_rotXYZ:
- alfa=2*M_PI*rando(&seed);
- beta=2*M_PI*rando(&seed);
- gamma=2*M_PI*rando(&seed);
- break;
- case en_rotZ:
- alfa=beta=0.;
- gamma=2*M_PI*rando(&seed);
- break;
- case en_rotNone:
- alfa=beta=gamma=0.;
- break;
+ gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
}
- if (enum_rot==en_rotXYZ || (enum_rot==en_rotZ))
+ /* allocate memory for atom coordinates of insert molecules */
+ snew(x_insrt, atoms_insrt.nr);
+ snew(r_insrt, atoms_insrt.nr);
+ snew(atoms_insrt.resinfo, atoms_insrt.nr);
+ snew(atoms_insrt.atomname, atoms_insrt.nr);
+ snew(atoms_insrt.atom, atoms_insrt.nr);
+ atoms_insrt.pdbinfo = NULL;
+ snew(x_n, atoms_insrt.nr);
+ snew(title_insrt, STRLEN);
+
+ /* read residue number, residue names, atomnames, coordinates etc. */
+ fprintf(stderr, "Reading molecule configuration \n");
+ read_stx_conf(mol_insrt, title_insrt, &atoms_insrt, x_insrt, NULL,
+ &ePBC_insrt, box_insrt);
+ fprintf(stderr, "%s\nContaining %d atoms in %d residue\n",
+ title_insrt, atoms_insrt.nr, atoms_insrt.nres);
+ srenew(atoms_insrt.resinfo, atoms_insrt.nres);
+
+ /* initialise van der waals arrays of insert molecules */
+ mk_vdw(&atoms_insrt, r_insrt, aps, r_distance);
+
+ /* With -ip, take nmol_insrt from file posfn */
+ if (posfn != NULL)
{
- rotate_conf(atoms_insrt.nr,x_n,NULL,alfa,beta,gamma);
+ nmol_insrt = read_xvg(posfn, &rpos, &ncol);
+ if (ncol != 3)
+ {
+ gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n", ncol, posfn);
+ }
+ fprintf(stderr, "Read %d positions from file %s\n\n", nmol_insrt, posfn);
}
- if (posfn==NULL)
+
+ srenew(atoms->resinfo, (atoms->nres+nmol_insrt*atoms_insrt.nres));
+ srenew(atoms->atomname, (atoms->nr+atoms_insrt.nr*nmol_insrt));
+ srenew(atoms->atom, (atoms->nr+atoms_insrt.nr*nmol_insrt));
+ srenew(*x, (atoms->nr+atoms_insrt.nr*nmol_insrt));
+ srenew(*r, (atoms->nr+atoms_insrt.nr*nmol_insrt));
+
+ trial = mol = 0;
+ while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
{
- /* insert at random positions */
- offset_x[XX]=box[XX][XX]*rando(&seed);
- offset_x[YY]=box[YY][YY]*rando(&seed);
- offset_x[ZZ]=box[ZZ][ZZ]*rando(&seed);
- gen_box(0,atoms_insrt.nr,x_n,box_insrt,offset_x,TRUE);
- if (!in_box(&pbc,x_n[0]) || !in_box(&pbc,x_n[atoms_insrt.nr-1]))
+ fprintf(stderr, "\rTry %d", trial++);
+ for (i = 0; (i < atoms_insrt.nr); i++)
+ {
+ copy_rvec(x_insrt[i], x_n[i]);
+ }
+ switch (enum_rot)
+ {
+ case en_rotXYZ:
+ alfa = 2*M_PI *rando(&seed);
+ beta = 2*M_PI *rando(&seed);
+ gamma = 2*M_PI *rando(&seed);
+ break;
+ case en_rotZ:
+ alfa = beta = 0.;
+ gamma = 2*M_PI*rando(&seed);
+ break;
+ case en_rotNone:
+ alfa = beta = gamma = 0.;
+ break;
+ }
+ if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
+ {
+ rotate_conf(atoms_insrt.nr, x_n, NULL, alfa, beta, gamma);
+ }
+ if (posfn == NULL)
+ {
+ /* insert at random positions */
+ offset_x[XX] = box[XX][XX]*rando(&seed);
+ offset_x[YY] = box[YY][YY]*rando(&seed);
+ offset_x[ZZ] = box[ZZ][ZZ]*rando(&seed);
+ gen_box(0, atoms_insrt.nr, x_n, box_insrt, offset_x, TRUE);
+ if (!in_box(&pbc, x_n[0]) || !in_box(&pbc, x_n[atoms_insrt.nr-1]))
+ {
+ continue;
+ }
+ }
+ else
+ {
+ /* Insert at positions taken from option -ip file */
+ offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2*rando(&seed)-1);
+ offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2*rando(&seed)-1);
+ offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2*rando(&seed)-1);
+ for (i = 0; i < atoms_insrt.nr; i++)
+ {
+ rvec_inc(x_n[i], offset_x);
+ }
+ }
+
+ onr = atoms->nr;
+
+ /* This is a (maybe) slow workaround to avoid too many calls of add_conf, which
+ * leaks memory (status May 2012). If the momory leaks in add_conf() are fixed,
+ * this check could be removed. Note, however, that allPairsDistOk is probably
+ * even faster than add_conf() when inserting a small molecule into a moderately
+ * small system.
+ */
+ if (bCheckAllPairDist && !allPairsDistOk(atoms, *x, *r, ePBC, box, &atoms_insrt, x_n, r_insrt))
+ {
continue;
+ }
+
+ add_conf(atoms, x, NULL, r, FALSE, ePBC, box, TRUE,
+ &atoms_insrt, x_n, NULL, r_insrt, FALSE, rshell, 0, oenv);
+
+ if (atoms->nr == (atoms_insrt.nr+onr))
+ {
+ mol++;
+ fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
+ }
}
- else
+ srenew(atoms->resinfo, atoms->nres);
+ srenew(atoms->atomname, atoms->nr);
+ srenew(atoms->atom, atoms->nr);
+ srenew(*x, atoms->nr);
+ srenew(*r, atoms->nr);
+
+ fprintf(stderr, "\n");
+ /* print number of molecules added */
+ fprintf(stderr, "Added %d molecules (out of %d requested) of %s\n",
+ mol, nmol_insrt, *atoms_insrt.resinfo[0].name);
+
+ return title_insrt;
+}
+
+static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **r,
+ int ePBC, matrix box,
+ gmx_atomprop_t aps, real r_distance, int *atoms_added,
+ int *residues_added, real rshell, int max_sol,
+ const output_env_t oenv)
+{
+ int i, nmol;
+ ivec n_box;
+ char filename[STRLEN];
+ char title_solvt[STRLEN];
+ t_atoms *atoms_solvt;
+ rvec *x_solvt, *v_solvt = NULL;
+ real *r_solvt;
+ int ePBC_solvt;
+ matrix box_solvt;
+ int onr, onres;
+ char *lfn;
+
+ lfn = gmxlibfn(fn);
+ strncpy(filename, lfn, STRLEN);
+ sfree(lfn);
+ snew(atoms_solvt, 1);
+ get_stx_coordnum(filename, &(atoms_solvt->nr));
+ if (atoms_solvt->nr == 0)
+ {
+ gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
+ }
+ snew(x_solvt, atoms_solvt->nr);
+ if (v)
{
- /* Insert at positions taken from option -ip file */
- offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2*rando(&seed)-1);
- offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2*rando(&seed)-1);
- offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2*rando(&seed)-1);
- for (i=0; i<atoms_insrt.nr; i++)
- rvec_inc(x_n[i], offset_x);
+ snew(v_solvt, atoms_solvt->nr);
}
+ snew(r_solvt, atoms_solvt->nr);
+ snew(atoms_solvt->resinfo, atoms_solvt->nr);
+ snew(atoms_solvt->atomname, atoms_solvt->nr);
+ snew(atoms_solvt->atom, atoms_solvt->nr);
+ atoms_solvt->pdbinfo = NULL;
+ fprintf(stderr, "Reading solvent configuration%s\n",
+ v_solvt ? " and velocities" : "");
+ read_stx_conf(filename, title_solvt, atoms_solvt, x_solvt, v_solvt,
+ &ePBC_solvt, box_solvt);
+ fprintf(stderr, "\"%s\"\n", title_solvt);
+ fprintf(stderr, "solvent configuration contains %d atoms in %d residues\n",
+ atoms_solvt->nr, atoms_solvt->nres);
+ fprintf(stderr, "\n");
- onr=atoms->nr;
+ /* apply pbc for solvent configuration for whole molecules */
+ rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
- /* This is a (maybe) slow workaround to avoid too many calls of add_conf, which
- * leaks memory (status May 2012). If the momory leaks in add_conf() are fixed,
- * this check could be removed. Note, however, that allPairsDistOk is probably
- * even faster than add_conf() when inserting a small molecule into a moderately
- * small system.
- */
- if (bCheckAllPairDist && !allPairsDistOk(atoms,*x,*r,ePBC,box,&atoms_insrt,x_n,r_insrt))
- continue;
+ /* initialise van der waals arrays of solvent configuration */
+ mk_vdw(atoms_solvt, r_solvt, aps, r_distance);
- add_conf(atoms,x,NULL,r,FALSE,ePBC,box,TRUE,
- &atoms_insrt,x_n,NULL,r_insrt,FALSE,rshell,0,oenv);
+ /* calculate the box multiplication factors n_box[0...DIM] */
+ nmol = 1;
+ for (i = 0; (i < DIM); i++)
+ {
+ n_box[i] = 1;
+ while (n_box[i]*box_solvt[i][i] < box[i][i])
+ {
+ n_box[i]++;
+ }
+ nmol *= n_box[i];
+ }
+ fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
+ n_box[XX], n_box[YY], n_box[ZZ]);
- if (atoms->nr==(atoms_insrt.nr+onr)) {
- mol++;
- fprintf(stderr," success (now %d atoms)!\n",atoms->nr);
+ /* realloc atoms_solvt for the new solvent configuration */
+ srenew(atoms_solvt->resinfo, atoms_solvt->nres*nmol);
+ srenew(atoms_solvt->atomname, atoms_solvt->nr*nmol);
+ srenew(atoms_solvt->atom, atoms_solvt->nr*nmol);
+ srenew(x_solvt, atoms_solvt->nr*nmol);
+ if (v_solvt)
+ {
+ srenew(v_solvt, atoms_solvt->nr*nmol);
}
- }
- srenew(atoms->resinfo, atoms->nres);
- srenew(atoms->atomname, atoms->nr);
- srenew(atoms->atom, atoms->nr);
- srenew(*x, atoms->nr);
- srenew(*r, atoms->nr);
-
- fprintf(stderr,"\n");
- /* print number of molecules added */
- fprintf(stderr,"Added %d molecules (out of %d requested) of %s\n",
- mol,nmol_insrt,*atoms_insrt.resinfo[0].name);
-
- return title_insrt;
-}
+ srenew(r_solvt, atoms_solvt->nr*nmol);
-static void add_solv(const char *fn,t_atoms *atoms,rvec **x,rvec **v,real **r,
- int ePBC,matrix box,
- gmx_atomprop_t aps,real r_distance,int *atoms_added,
- int *residues_added,real rshell,int max_sol,
- const output_env_t oenv)
-{
- int i,nmol;
- ivec n_box;
- char filename[STRLEN];
- char title_solvt[STRLEN];
- t_atoms *atoms_solvt;
- rvec *x_solvt,*v_solvt=NULL;
- real *r_solvt;
- int ePBC_solvt;
- matrix box_solvt;
- int onr,onres;
- char *lfn;
-
- lfn = gmxlibfn(fn);
- strncpy(filename,lfn,STRLEN);
- sfree(lfn);
- snew(atoms_solvt,1);
- get_stx_coordnum(filename,&(atoms_solvt->nr));
- if (atoms_solvt->nr == 0)
- gmx_fatal(FARGS,"No solvent in %s, please check your input\n",filename);
- snew(x_solvt,atoms_solvt->nr);
- if (v) snew(v_solvt,atoms_solvt->nr);
- snew(r_solvt,atoms_solvt->nr);
- snew(atoms_solvt->resinfo,atoms_solvt->nr);
- snew(atoms_solvt->atomname,atoms_solvt->nr);
- snew(atoms_solvt->atom,atoms_solvt->nr);
- atoms_solvt->pdbinfo = NULL;
- fprintf(stderr,"Reading solvent configuration%s\n",
- v_solvt?" and velocities":"");
- read_stx_conf(filename,title_solvt,atoms_solvt,x_solvt,v_solvt,
- &ePBC_solvt,box_solvt);
- fprintf(stderr,"\"%s\"\n",title_solvt);
- fprintf(stderr,"solvent configuration contains %d atoms in %d residues\n",
- atoms_solvt->nr,atoms_solvt->nres);
- fprintf(stderr,"\n");
-
- /* apply pbc for solvent configuration for whole molecules */
- rm_res_pbc(atoms_solvt,x_solvt,box_solvt);
-
- /* initialise van der waals arrays of solvent configuration */
- mk_vdw(atoms_solvt,r_solvt,aps,r_distance);
-
- /* calculate the box multiplication factors n_box[0...DIM] */
- nmol=1;
- for (i=0; (i < DIM);i++) {
- n_box[i] = 1;
- while (n_box[i]*box_solvt[i][i] < box[i][i])
- n_box[i]++;
- nmol*=n_box[i];
- }
- fprintf(stderr,"Will generate new solvent configuration of %dx%dx%d boxes\n",
- n_box[XX],n_box[YY],n_box[ZZ]);
-
- /* realloc atoms_solvt for the new solvent configuration */
- srenew(atoms_solvt->resinfo,atoms_solvt->nres*nmol);
- srenew(atoms_solvt->atomname,atoms_solvt->nr*nmol);
- srenew(atoms_solvt->atom,atoms_solvt->nr*nmol);
- srenew(x_solvt,atoms_solvt->nr*nmol);
- if (v_solvt) srenew(v_solvt,atoms_solvt->nr*nmol);
- srenew(r_solvt,atoms_solvt->nr*nmol);
-
- /* generate a new solvent configuration */
- genconf(atoms_solvt,x_solvt,v_solvt,r_solvt,box_solvt,n_box);
+ /* generate a new solvent configuration */
+ genconf(atoms_solvt, x_solvt, v_solvt, r_solvt, box_solvt, n_box);
#ifdef DEBUG
- print_stat(x_solvt,atoms_solvt->nr,box_solvt);
+ print_stat(x_solvt, atoms_solvt->nr, box_solvt);
#endif
-
+
#ifdef DEBUG
- print_stat(x_solvt,atoms_solvt->nr,box_solvt);
+ print_stat(x_solvt, atoms_solvt->nr, box_solvt);
#endif
- /* Sort the solvent mixture, not the protein... */
- sort_molecule(&atoms_solvt,x_solvt,v_solvt,r_solvt);
-
- /* add the two configurations */
- onr=atoms->nr;
- onres=atoms->nres;
- add_conf(atoms,x,v,r,TRUE,ePBC,box,FALSE,
- atoms_solvt,x_solvt,v_solvt,r_solvt,TRUE,rshell,max_sol,oenv);
- *atoms_added=atoms->nr-onr;
- *residues_added=atoms->nres-onres;
-
- sfree(x_solvt);
- sfree(r_solvt);
-
- fprintf(stderr,"Generated solvent containing %d atoms in %d residues\n",
- *atoms_added,*residues_added);
+ /* Sort the solvent mixture, not the protein... */
+ sort_molecule(&atoms_solvt, x_solvt, v_solvt, r_solvt);
+
+ /* add the two configurations */
+ onr = atoms->nr;
+ onres = atoms->nres;
+ add_conf(atoms, x, v, r, TRUE, ePBC, box, FALSE,
+ atoms_solvt, x_solvt, v_solvt, r_solvt, TRUE, rshell, max_sol, oenv);
+ *atoms_added = atoms->nr-onr;
+ *residues_added = atoms->nres-onres;
+
+ sfree(x_solvt);
+ sfree(r_solvt);
+
+ fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
+ *atoms_added, *residues_added);
}
-static char *read_prot(const char *confin,t_atoms *atoms,rvec **x,rvec **v,
- real **r, int *ePBC,matrix box,gmx_atomprop_t aps,
+static char *read_prot(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
+ real **r, int *ePBC, matrix box, gmx_atomprop_t aps,
real r_distance)
{
- char *title;
- int natoms;
-
- snew(title,STRLEN);
- get_stx_coordnum(confin,&natoms);
-
- /* allocate memory for atom coordinates of configuration 1 */
- snew(*x,natoms);
- if (v) snew(*v,natoms);
- snew(*r,natoms);
- init_t_atoms(atoms,natoms,FALSE);
-
- /* read residue number, residue names, atomnames, coordinates etc. */
- fprintf(stderr,"Reading solute configuration%s\n",v?" and velocities":"");
- read_stx_conf(confin,title,atoms,*x,v?*v:NULL,ePBC,box);
- fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
- title,atoms->nr,atoms->nres);
-
- /* initialise van der waals arrays of configuration 1 */
- mk_vdw(atoms,*r,aps,r_distance);
-
- return title;
+ char *title;
+ int natoms;
+
+ snew(title, STRLEN);
+ get_stx_coordnum(confin, &natoms);
+
+ /* allocate memory for atom coordinates of configuration 1 */
+ snew(*x, natoms);
+ if (v)
+ {
+ snew(*v, natoms);
+ }
+ snew(*r, natoms);
+ init_t_atoms(atoms, natoms, FALSE);
+
+ /* read residue number, residue names, atomnames, coordinates etc. */
+ fprintf(stderr, "Reading solute configuration%s\n", v ? " and velocities" : "");
+ read_stx_conf(confin, title, atoms, *x, v ? *v : NULL, ePBC, box);
+ fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
+ title, atoms->nr, atoms->nres);
+
+ /* initialise van der waals arrays of configuration 1 */
+ mk_vdw(atoms, *r, aps, r_distance);
+
+ return title;
}
-static void update_top(t_atoms *atoms,matrix box,int NFILE,t_filenm fnm[],
- gmx_atomprop_t aps)
+static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
+ gmx_atomprop_t aps)
{
#define TEMP_FILENM "temp.top"
- FILE *fpin,*fpout;
- char buf[STRLEN],buf2[STRLEN],*temp;
- const char *topinout;
- int line;
- gmx_bool bSystem,bMolecules,bSkip;
- int i,nsol=0;
- double mtot;
- real vol,mm;
-
- for(i=0; (i<atoms->nres); i++) {
- /* calculate number of SOLvent molecules */
- if ( (strcmp(*atoms->resinfo[i].name,"SOL")==0) ||
- (strcmp(*atoms->resinfo[i].name,"WAT")==0) ||
- (strcmp(*atoms->resinfo[i].name,"HOH")==0) )
- nsol++;
- }
- mtot = 0;
- for(i=0; (i<atoms->nr); i++) {
- gmx_atomprop_query(aps,epropMass,
- *atoms->resinfo[atoms->atom[i].resind].name,
- *atoms->atomname[i],&mm);
- mtot += mm;
- }
-
- vol=det(box);
-
- fprintf(stderr,"Volume : %10g (nm^3)\n",vol);
- fprintf(stderr,"Density : %10g (g/l)\n",
- (mtot*1e24)/(AVOGADRO*vol));
- fprintf(stderr,"Number of SOL molecules: %5d \n\n",nsol);
-
- /* open topology file and append sol molecules */
- topinout = ftp2fn(efTOP,NFILE,fnm);
- if ( ftp2bSet(efTOP,NFILE,fnm) ) {
- fprintf(stderr,"Processing topology\n");
- fpin = ffopen(topinout,"r");
- fpout= ffopen(TEMP_FILENM,"w");
- line=0;
- bSystem = bMolecules = FALSE;
- while (fgets(buf, STRLEN, fpin)) {
- bSkip=FALSE;
- line++;
- strcpy(buf2,buf);
- if ((temp=strchr(buf2,'\n')) != NULL)
- temp[0]='\0';
- ltrim(buf2);
- if (buf2[0]=='[') {
- buf2[0]=' ';
- if ((temp=strchr(buf2,'\n')) != NULL)
- temp[0]='\0';
- rtrim(buf2);
- if (buf2[strlen(buf2)-1]==']') {
- buf2[strlen(buf2)-1]='\0';
- ltrim(buf2);
- rtrim(buf2);
- bSystem=(gmx_strcasecmp(buf2,"system")==0);
- bMolecules=(gmx_strcasecmp(buf2,"molecules")==0);
- }
- } else if (bSystem && nsol && (buf[0]!=';') ) {
- /* if sol present, append "in water" to system name */
- rtrim(buf2);
- if (buf2[0] && (!strstr(buf2," water")) ) {
- sprintf(buf,"%s in water\n",buf2);
- bSystem = FALSE;
- }
- } else if (bMolecules) {
- /* check if this is a line with solvent molecules */
- sscanf(buf,"%4095s",buf2);
- if (strcmp(buf2,"SOL")==0) {
- sscanf(buf,"%*4095s %20d",&i);
- nsol-=i;
- if (nsol<0) {
- bSkip=TRUE;
- nsol+=i;
- }
- }
- }
- if (bSkip) {
- if ((temp=strchr(buf,'\n')) != NULL)
- temp[0]='\0';
- fprintf(stdout,"Removing line #%d '%s' from topology file (%s)\n",
- line,buf,topinout);
- } else
- fprintf(fpout,"%s",buf);
+ FILE *fpin, *fpout;
+ char buf[STRLEN], buf2[STRLEN], *temp;
+ const char *topinout;
+ int line;
+ gmx_bool bSystem, bMolecules, bSkip;
+ int i, nsol = 0;
+ double mtot;
+ real vol, mm;
+
+ for (i = 0; (i < atoms->nres); i++)
+ {
+ /* calculate number of SOLvent molecules */
+ if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
+ (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
+ (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
+ {
+ nsol++;
+ }
+ }
+ mtot = 0;
+ for (i = 0; (i < atoms->nr); i++)
+ {
+ gmx_atomprop_query(aps, epropMass,
+ *atoms->resinfo[atoms->atom[i].resind].name,
+ *atoms->atomname[i], &mm);
+ mtot += mm;
}
- ffclose(fpin);
- if ( nsol ) {
- fprintf(stdout,"Adding line for %d solvent molecules to "
- "topology file (%s)\n",nsol,topinout);
- fprintf(fpout,"%-15s %5d\n","SOL",nsol);
+
+ vol = det(box);
+
+ fprintf(stderr, "Volume : %10g (nm^3)\n", vol);
+ fprintf(stderr, "Density : %10g (g/l)\n",
+ (mtot*1e24)/(AVOGADRO*vol));
+ fprintf(stderr, "Number of SOL molecules: %5d \n\n", nsol);
+
+ /* open topology file and append sol molecules */
+ topinout = ftp2fn(efTOP, NFILE, fnm);
+ if (ftp2bSet(efTOP, NFILE, fnm) )
+ {
+ fprintf(stderr, "Processing topology\n");
+ fpin = ffopen(topinout, "r");
+ fpout = ffopen(TEMP_FILENM, "w");
+ line = 0;
+ bSystem = bMolecules = FALSE;
+ while (fgets(buf, STRLEN, fpin))
+ {
+ bSkip = FALSE;
+ line++;
+ strcpy(buf2, buf);
+ if ((temp = strchr(buf2, '\n')) != NULL)
+ {
+ temp[0] = '\0';
+ }
+ ltrim(buf2);
+ if (buf2[0] == '[')
+ {
+ buf2[0] = ' ';
+ if ((temp = strchr(buf2, '\n')) != NULL)
+ {
+ temp[0] = '\0';
+ }
+ rtrim(buf2);
+ if (buf2[strlen(buf2)-1] == ']')
+ {
+ buf2[strlen(buf2)-1] = '\0';
+ ltrim(buf2);
+ rtrim(buf2);
+ bSystem = (gmx_strcasecmp(buf2, "system") == 0);
+ bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
+ }
+ }
+ else if (bSystem && nsol && (buf[0] != ';') )
+ {
+ /* if sol present, append "in water" to system name */
+ rtrim(buf2);
+ if (buf2[0] && (!strstr(buf2, " water")) )
+ {
+ sprintf(buf, "%s in water\n", buf2);
+ bSystem = FALSE;
+ }
+ }
+ else if (bMolecules)
+ {
+ /* check if this is a line with solvent molecules */
+ sscanf(buf, "%4095s", buf2);
+ if (strcmp(buf2, "SOL") == 0)
+ {
+ sscanf(buf, "%*4095s %20d", &i);
+ nsol -= i;
+ if (nsol < 0)
+ {
+ bSkip = TRUE;
+ nsol += i;
+ }
+ }
+ }
+ if (bSkip)
+ {
+ if ((temp = strchr(buf, '\n')) != NULL)
+ {
+ temp[0] = '\0';
+ }
+ fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
+ line, buf, topinout);
+ }
+ else
+ {
+ fprintf(fpout, "%s", buf);
+ }
+ }
+ ffclose(fpin);
+ if (nsol)
+ {
+ fprintf(stdout, "Adding line for %d solvent molecules to "
+ "topology file (%s)\n", nsol, topinout);
+ fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
+ }
+ ffclose(fpout);
+ /* use ffopen to generate backup of topinout */
+ fpout = ffopen(topinout, "w");
+ ffclose(fpout);
+ rename(TEMP_FILENM, topinout);
}
- ffclose(fpout);
- /* use ffopen to generate backup of topinout */
- fpout=ffopen(topinout,"w");
- ffclose(fpout);
- rename(TEMP_FILENM,topinout);
- }
#undef TEMP_FILENM
}
-int gmx_genbox(int argc,char *argv[])
+int gmx_genbox(int argc, char *argv[])
{
- const char *desc[] = {
- "[TT]genbox[tt] can do one of 4 things:[PAR]",
-
- "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt]. Or specify [TT]-cs[tt] and",
- "[TT]-cp[tt] with a structure file with a box, but without atoms.[PAR]",
-
- "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
- "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
- "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
- "unless [TT]-box[tt] is set.",
- "If you want the solute to be centered in the box,",
- "the program [TT]editconf[tt] has sophisticated options",
- "to change the box dimensions and center the solute.",
- "Solvent molecules are removed from the box where the ",
- "distance between any atom of the solute molecule(s) and any atom of ",
- "the solvent molecule is less than the sum of the van der Waals radii of ",
- "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
- "read by the program, and atoms not in the database are ",
- "assigned a default distance [TT]-vdwd[tt].",
- "Note that this option will also influence the distances between",
- "solvent molecules if they contain atoms that are not in the database.",
- "[PAR]",
-
- "3) Insert a number ([TT]-nmol[tt]) of extra molecules ([TT]-ci[tt]) ",
- "at random positions.",
- "The program iterates until [TT]nmol[tt] molecules",
- "have been inserted in the box. To test whether an insertion is ",
- "successful the same van der Waals criterium is used as for removal of ",
- "solvent molecules. When no appropriately-sized ",
- "holes (holes that can hold an extra molecule) are available, the ",
- "program tries for [TT]-nmol[tt] * [TT]-try[tt] times before giving up. ",
- "Increase [TT]-try[tt] if you have several small holes to fill.",
- "Option [TT]-rot[tt] defines if the molecules are randomly oriented.",
- "[PAR]",
-
- "4) Insert a number of molecules ([TT]-ci[tt]) at positions defined in",
- "positions.dat ([TT]-ip[tt]). positions.dat should have 3 columns (x/y/z),",
- "that give the displacements compared to the input molecule position ([TT]-ci[tt]).",
- "Hence, if positions.dat should contain the absolut positions, the molecule ",
- "must be centered to 0/0/0 before using genbox (use, e.g., editconf -center).",
- "Comments in positions.dat starting with # are ignored. Option [TT]-dr[tt]",
- "defines the maximally allowed displacements during insertial trials.",
- "[TT]-try[tt] and [TT]-rot[tt] work as in mode (3) (see above)",
- "[PAR]",
-
- "If you need to do more than one of the above operations, it can be",
- "best to call [TT]genbox[tt] separately for each operation, so that",
- "you are sure of the order in which the operations occur.[PAR]",
-
- "The default solvent is Simple Point Charge water (SPC), with coordinates ",
- "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
- "for other 3-site water models, since a short equibilibration will remove",
- "the small differences between the models.",
- "Other solvents are also supported, as well as mixed solvents. The",
- "only restriction to solvent types is that a solvent molecule consists",
- "of exactly one residue. The residue information in the coordinate",
- "files is used, and should therefore be more or less consistent.",
- "In practice this means that two subsequent solvent molecules in the ",
- "solvent coordinate file should have different residue number.",
- "The box of solute is built by stacking the coordinates read from",
- "the coordinate file. This means that these coordinates should be ",
- "equlibrated in periodic boundary conditions to ensure a good",
- "alignment of molecules on the stacking interfaces.",
- "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
- "solvent molecules and leaves out the rest that would have fitted",
- "into the box. This can create a void that can cause problems later.",
- "Choose your volume wisely.[PAR]",
-
- "The program can optionally rotate the solute molecule to align the",
- "longest molecule axis along a box edge. This way the amount of solvent",
- "molecules necessary is reduced.",
- "It should be kept in mind that this only works for",
- "short simulations, as e.g. an alpha-helical peptide in solution can ",
- "rotate over 90 degrees, within 500 ps. In general it is therefore ",
- "better to make a more or less cubic box.[PAR]",
-
- "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
- "the specified thickness (nm) around the solute. Hint: it is a good",
- "idea to put the protein in the center of a box first (using [TT]editconf[tt]).",
- "[PAR]",
-
- "Finally, [TT]genbox[tt] will optionally remove lines from your topology file in ",
- "which a number of solvent molecules is already added, and adds a ",
- "line with the total number of solvent molecules in your coordinate file."
- };
-
- const char *bugs[] = {
- "Molecules must be whole in the initial configurations.",
- "Many repeated neighbor searchings with -ci blows up the allocated memory. "
- "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
- };
-
- /* parameter data */
- gmx_bool bSol,bProt,bBox;
- const char *conf_prot,*confout;
- int bInsert;
- real *r;
- char *title_ins;
- gmx_atomprop_t aps;
-
- /* protein configuration data */
- char *title=NULL;
- t_atoms atoms;
- rvec *x,*v=NULL;
- int ePBC=-1;
- matrix box;
-
- /* other data types */
- int atoms_added,residues_added;
-
- t_filenm fnm[] = {
- { efSTX, "-cp", "protein", ffOPTRD },
- { efSTX, "-cs", "spc216", ffLIBOPTRD},
- { efSTX, "-ci", "insert", ffOPTRD},
- { efDAT, "-ip", "positions", ffOPTRD},
- { efSTO, NULL, NULL, ffWRITE},
- { efTOP, NULL, NULL, ffOPTRW},
- };
+ const char *desc[] = {
+ "[TT]genbox[tt] can do one of 4 things:[PAR]",
+
+ "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt]. Or specify [TT]-cs[tt] and",
+ "[TT]-cp[tt] with a structure file with a box, but without atoms.[PAR]",
+
+ "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
+ "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
+ "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
+ "unless [TT]-box[tt] is set.",
+ "If you want the solute to be centered in the box,",
+ "the program [TT]editconf[tt] has sophisticated options",
+ "to change the box dimensions and center the solute.",
+ "Solvent molecules are removed from the box where the ",
+ "distance between any atom of the solute molecule(s) and any atom of ",
+ "the solvent molecule is less than the sum of the van der Waals radii of ",
+ "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
+ "read by the program, and atoms not in the database are ",
+ "assigned a default distance [TT]-vdwd[tt].",
+ "Note that this option will also influence the distances between",
+ "solvent molecules if they contain atoms that are not in the database.",
+ "[PAR]",
+
+ "3) Insert a number ([TT]-nmol[tt]) of extra molecules ([TT]-ci[tt]) ",
+ "at random positions.",
+ "The program iterates until [TT]nmol[tt] molecules",
+ "have been inserted in the box. To test whether an insertion is ",
+ "successful the same van der Waals criterium is used as for removal of ",
+ "solvent molecules. When no appropriately-sized ",
+ "holes (holes that can hold an extra molecule) are available, the ",
+ "program tries for [TT]-nmol[tt] * [TT]-try[tt] times before giving up. ",
+ "Increase [TT]-try[tt] if you have several small holes to fill.",
+ "Option [TT]-rot[tt] defines if the molecules are randomly oriented.",
+ "[PAR]",
+
+ "4) Insert a number of molecules ([TT]-ci[tt]) at positions defined in",
+ "positions.dat ([TT]-ip[tt]). positions.dat should have 3 columns (x/y/z),",
+ "that give the displacements compared to the input molecule position ([TT]-ci[tt]).",
+ "Hence, if positions.dat should contain the absolut positions, the molecule ",
+ "must be centered to 0/0/0 before using genbox (use, e.g., editconf -center).",
+ "Comments in positions.dat starting with # are ignored. Option [TT]-dr[tt]",
+ "defines the maximally allowed displacements during insertial trials.",
+ "[TT]-try[tt] and [TT]-rot[tt] work as in mode (3) (see above)",
+ "[PAR]",
+
+ "If you need to do more than one of the above operations, it can be",
+ "best to call [TT]genbox[tt] separately for each operation, so that",
+ "you are sure of the order in which the operations occur.[PAR]",
+
+ "The default solvent is Simple Point Charge water (SPC), with coordinates ",
+ "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
+ "for other 3-site water models, since a short equibilibration will remove",
+ "the small differences between the models.",
+ "Other solvents are also supported, as well as mixed solvents. The",
+ "only restriction to solvent types is that a solvent molecule consists",
+ "of exactly one residue. The residue information in the coordinate",
+ "files is used, and should therefore be more or less consistent.",
+ "In practice this means that two subsequent solvent molecules in the ",
+ "solvent coordinate file should have different residue number.",
+ "The box of solute is built by stacking the coordinates read from",
+ "the coordinate file. This means that these coordinates should be ",
+ "equlibrated in periodic boundary conditions to ensure a good",
+ "alignment of molecules on the stacking interfaces.",
+ "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
+ "solvent molecules and leaves out the rest that would have fitted",
+ "into the box. This can create a void that can cause problems later.",
+ "Choose your volume wisely.[PAR]",
+
+ "The program can optionally rotate the solute molecule to align the",
+ "longest molecule axis along a box edge. This way the amount of solvent",
+ "molecules necessary is reduced.",
+ "It should be kept in mind that this only works for",
+ "short simulations, as e.g. an alpha-helical peptide in solution can ",
+ "rotate over 90 degrees, within 500 ps. In general it is therefore ",
+ "better to make a more or less cubic box.[PAR]",
+
+ "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
+ "the specified thickness (nm) around the solute. Hint: it is a good",
+ "idea to put the protein in the center of a box first (using [TT]editconf[tt]).",
+ "[PAR]",
+
+ "Finally, [TT]genbox[tt] will optionally remove lines from your topology file in ",
+ "which a number of solvent molecules is already added, and adds a ",
+ "line with the total number of solvent molecules in your coordinate file."
+ };
+
+ const char *bugs[] = {
+ "Molecules must be whole in the initial configurations.",
+ "Many repeated neighbor searchings with -ci blows up the allocated memory. "
+ "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
+ };
+
+ /* parameter data */
+ gmx_bool bSol, bProt, bBox;
+ const char *conf_prot, *confout;
+ int bInsert;
+ real *r;
+ char *title_ins;
+ gmx_atomprop_t aps;
+
+ /* protein configuration data */
+ char *title = NULL;
+ t_atoms atoms;
+ rvec *x, *v = NULL;
+ int ePBC = -1;
+ matrix box;
+
+ /* other data types */
+ int atoms_added, residues_added;
+
+ t_filenm fnm[] = {
+ { efSTX, "-cp", "protein", ffOPTRD },
+ { efSTX, "-cs", "spc216", ffLIBOPTRD},
+ { efSTX, "-ci", "insert", ffOPTRD},
+ { efDAT, "-ip", "positions", ffOPTRD},
+ { efSTO, NULL, NULL, ffWRITE},
+ { efTOP, NULL, NULL, ffOPTRW},
+ };
#define NFILE asize(fnm)
-
- static int nmol_ins=0,nmol_try=10,seed=1997,enum_rot;
- static real r_distance=0.105,r_shell=0;
- static rvec new_box={0.0,0.0,0.0}, deltaR={0.0,0.0,0.0};
- static gmx_bool bReadV=FALSE,bCheckAllPairDist=FALSE;
- static int max_sol = 0;
- output_env_t oenv;
- const char *enum_rot_string[]={NULL,"xyz","z","none",NULL};
- t_pargs pa[] = {
- { "-box", FALSE, etRVEC, {new_box},
- "Box size" },
- { "-nmol", FALSE, etINT , {&nmol_ins},
- "Number of extra molecules to insert" },
- { "-try", FALSE, etINT , {&nmol_try},
- "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
- { "-seed", FALSE, etINT , {&seed},
- "Random generator seed"},
- { "-vdwd", FALSE, etREAL, {&r_distance},
- "Default van der Waals distance"},
- { "-shell", FALSE, etREAL, {&r_shell},
- "Thickness of optional water layer around solute" },
- { "-maxsol", FALSE, etINT, {&max_sol},
- "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
- { "-vel", FALSE, etBOOL, {&bReadV},
- "Keep velocities from input solute and solvent" },
- { "-dr", FALSE, etRVEC, {deltaR},
- "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
- { "-rot", FALSE, etENUM, {enum_rot_string},
- "rotate inserted molecules randomly" },
- { "-allpair", FALSE, etBOOL, {&bCheckAllPairDist},
- "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
- };
-
- parse_common_args(&argc,argv, PCA_BE_NICE,NFILE,fnm,asize(pa),pa,
- asize(desc),desc,asize(bugs),bugs,&oenv);
-
- bInsert = opt2bSet("-ci",NFILE,fnm) && ((nmol_ins > 0) || opt2bSet("-ip",NFILE,fnm));
- bSol = opt2bSet("-cs",NFILE,fnm);
- bProt = opt2bSet("-cp",NFILE,fnm);
- bBox = opt2parg_bSet("-box",asize(pa),pa);
- enum_rot = nenum(enum_rot_string);
-
- /* check input */
- if (bInsert && (nmol_ins<=0 && !opt2bSet("-ip",NFILE,fnm)))
- gmx_fatal(FARGS,"When specifying inserted molecules (-ci), "
- "-nmol must be larger than 0 or positions must be given with -ip");
- if (!bProt && !bBox)
- gmx_fatal(FARGS,"When no solute (-cp) is specified, "
- "a box size (-box) must be specified");
-
- aps = gmx_atomprop_init();
-
- if (bProt) {
- /*generate a solute configuration */
- conf_prot = opt2fn("-cp",NFILE,fnm);
- title = read_prot(conf_prot,&atoms,&x,bReadV?&v:NULL,&r,&ePBC,box,
- aps,r_distance);
- if (bReadV && !v)
- fprintf(stderr,"Note: no velocities found\n");
- if (atoms.nr == 0) {
- fprintf(stderr,"Note: no atoms in %s\n",conf_prot);
- bProt = FALSE;
+
+ static int nmol_ins = 0, nmol_try = 10, seed = 1997, enum_rot;
+ static real r_distance = 0.105, r_shell = 0;
+ static rvec new_box = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
+ static gmx_bool bReadV = FALSE, bCheckAllPairDist = FALSE;
+ static int max_sol = 0;
+ output_env_t oenv;
+ const char *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
+ t_pargs pa[] = {
+ { "-box", FALSE, etRVEC, {new_box},
+ "Box size" },
+ { "-nmol", FALSE, etINT, {&nmol_ins},
+ "Number of extra molecules to insert" },
+ { "-try", FALSE, etINT, {&nmol_try},
+ "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
+ { "-seed", FALSE, etINT, {&seed},
+ "Random generator seed"},
+ { "-vdwd", FALSE, etREAL, {&r_distance},
+ "Default van der Waals distance"},
+ { "-shell", FALSE, etREAL, {&r_shell},
+ "Thickness of optional water layer around solute" },
+ { "-maxsol", FALSE, etINT, {&max_sol},
+ "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
+ { "-vel", FALSE, etBOOL, {&bReadV},
+ "Keep velocities from input solute and solvent" },
+ { "-dr", FALSE, etRVEC, {deltaR},
+ "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
+ { "-rot", FALSE, etENUM, {enum_rot_string},
+ "rotate inserted molecules randomly" },
+ { "-allpair", FALSE, etBOOL, {&bCheckAllPairDist},
+ "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
+ };
+
+ parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
+ asize(desc), desc, asize(bugs), bugs, &oenv);
+
+ bInsert = opt2bSet("-ci", NFILE, fnm) && ((nmol_ins > 0) || opt2bSet("-ip", NFILE, fnm));
+ bSol = opt2bSet("-cs", NFILE, fnm);
+ bProt = opt2bSet("-cp", NFILE, fnm);
+ bBox = opt2parg_bSet("-box", asize(pa), pa);
+ enum_rot = nenum(enum_rot_string);
+
+ /* check input */
+ if (bInsert && (nmol_ins <= 0 && !opt2bSet("-ip", NFILE, fnm)))
+ {
+ gmx_fatal(FARGS, "When specifying inserted molecules (-ci), "
+ "-nmol must be larger than 0 or positions must be given with -ip");
+ }
+ if (!bProt && !bBox)
+ {
+ gmx_fatal(FARGS, "When no solute (-cp) is specified, "
+ "a box size (-box) must be specified");
+ }
+
+ aps = gmx_atomprop_init();
+
+ if (bProt)
+ {
+ /*generate a solute configuration */
+ conf_prot = opt2fn("-cp", NFILE, fnm);
+ title = read_prot(conf_prot, &atoms, &x, bReadV ? &v : NULL, &r, &ePBC, box,
+ aps, r_distance);
+ if (bReadV && !v)
+ {
+ fprintf(stderr, "Note: no velocities found\n");
+ }
+ if (atoms.nr == 0)
+ {
+ fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
+ bProt = FALSE;
+ }
+ }
+ if (!bProt)
+ {
+ atoms.nr = 0;
+ atoms.nres = 0;
+ atoms.resinfo = NULL;
+ atoms.atomname = NULL;
+ atoms.atom = NULL;
+ atoms.pdbinfo = NULL;
+ x = NULL;
+ r = NULL;
+ }
+ if (bBox)
+ {
+ ePBC = epbcXYZ;
+ clear_mat(box);
+ box[XX][XX] = new_box[XX];
+ box[YY][YY] = new_box[YY];
+ box[ZZ][ZZ] = new_box[ZZ];
}
- }
- if (!bProt) {
- atoms.nr=0;
- atoms.nres=0;
- atoms.resinfo=NULL;
- atoms.atomname=NULL;
- atoms.atom=NULL;
- atoms.pdbinfo=NULL;
- x=NULL;
- r=NULL;
- }
- if (bBox) {
- ePBC = epbcXYZ;
- clear_mat(box);
- box[XX][XX]=new_box[XX];
- box[YY][YY]=new_box[YY];
- box[ZZ][ZZ]=new_box[ZZ];
- }
- if (det(box) == 0)
- gmx_fatal(FARGS,"Undefined solute box.\nCreate one with editconf "
- "or give explicit -box command line option");
-
- /* add nmol_ins molecules of atoms_ins
- in random orientation at random place */
- if (bInsert)
- title_ins = insert_mols(opt2fn("-ci",NFILE,fnm),nmol_ins,nmol_try,seed,
- &atoms,&x,&r,ePBC,box,aps,r_distance,r_shell,
- oenv,opt2fn_null("-ip",NFILE,fnm),deltaR,enum_rot,
- bCheckAllPairDist);
- else
- title_ins = strdup("Generated by genbox");
-
- /* add solvent */
- if (bSol)
- add_solv(opt2fn("-cs",NFILE,fnm),&atoms,&x,v?&v:NULL,&r,ePBC,box,
- aps,r_distance,&atoms_added,&residues_added,r_shell,max_sol,
- oenv);
-
- /* write new configuration 1 to file confout */
- confout = ftp2fn(efSTO,NFILE,fnm);
- fprintf(stderr,"Writing generated configuration to %s\n",confout);
- if (bProt) {
- write_sto_conf(confout,title,&atoms,x,v,ePBC,box);
- /* print box sizes and box type to stderr */
- fprintf(stderr,"%s\n",title);
- } else
- write_sto_conf(confout,title_ins,&atoms,x,v,ePBC,box);
-
- sfree(title_ins);
-
- /* print size of generated configuration */
- fprintf(stderr,"\nOutput configuration contains %d atoms in %d residues\n",
- atoms.nr,atoms.nres);
- update_top(&atoms,box,NFILE,fnm,aps);
-
- gmx_atomprop_destroy(aps);
-
- thanx(stderr);
-
- return 0;
+ if (det(box) == 0)
+ {
+ gmx_fatal(FARGS, "Undefined solute box.\nCreate one with editconf "
+ "or give explicit -box command line option");
+ }
+
+ /* add nmol_ins molecules of atoms_ins
+ in random orientation at random place */
+ if (bInsert)
+ {
+ title_ins = insert_mols(opt2fn("-ci", NFILE, fnm), nmol_ins, nmol_try, seed,
+ &atoms, &x, &r, ePBC, box, aps, r_distance, r_shell,
+ oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
+ bCheckAllPairDist);
+ }
+ else
+ {
+ title_ins = strdup("Generated by genbox");
+ }
+
+ /* add solvent */
+ if (bSol)
+ {
+ add_solv(opt2fn("-cs", NFILE, fnm), &atoms, &x, v ? &v : NULL, &r, ePBC, box,
+ aps, r_distance, &atoms_added, &residues_added, r_shell, max_sol,
+ oenv);
+ }
+
+ /* write new configuration 1 to file confout */
+ confout = ftp2fn(efSTO, NFILE, fnm);
+ fprintf(stderr, "Writing generated configuration to %s\n", confout);
+ if (bProt)
+ {
+ write_sto_conf(confout, title, &atoms, x, v, ePBC, box);
+ /* print box sizes and box type to stderr */
+ fprintf(stderr, "%s\n", title);
+ }
+ else
+ {
+ write_sto_conf(confout, title_ins, &atoms, x, v, ePBC, box);
+ }
+
+ sfree(title_ins);
+
+ /* print size of generated configuration */
+ fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
+ atoms.nr, atoms.nres);
+ update_top(&atoms, box, NFILE, fnm, aps);
+
+ gmx_atomprop_destroy(aps);
+
+ thanx(stderr);
+
+ return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff