* allPairsDistOk is probably even faster than the other code.
*/
static gmx_bool
-allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
+allPairsDistOk(t_atoms *atoms, rvec *x, real *exclusionDistances,
int ePBC, matrix box,
- t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
+ t_atoms *atoms_insrt, rvec *x_n, real *exclusionDistances_insrt)
{
int i, j;
rvec dx;
{
pbc_dx(&pbc, x[i], x_n[j], dx);
n2 = norm2(dx);
- r2 = sqr(r[i]+r_insrt[j]);
+ r2 = sqr(exclusionDistances[i]+exclusionDistances_insrt[j]);
if (n2 < r2)
{
return FALSE;
};
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,
+ t_atoms *atoms, rvec **x, real **exclusionDistances, int ePBC, matrix box,
gmx_atomprop_t aps,
- real r_distance, real r_scale, real rshell,
+ real defaultDistance, real scaleFactor, real rshell,
const output_env_t oenv,
const char* posfn, const rvec deltaR, int enum_rot,
gmx_bool bCheckAllPairDist)
static char *title_insrt;
t_atoms atoms_insrt;
rvec *x_insrt, *x_n;
- real *r_insrt;
+ real *exclusionDistances_insrt;
int ePBC_insrt;
matrix box_insrt;
int i, mol, onr, ncol;
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);
+ int natoms;
+ get_stx_coordnum(mol_insrt, &natoms);
+ if (natoms == 0)
+ {
+ gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
+ }
+ init_t_atoms(&atoms_insrt, natoms, FALSE);
}
/* 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);
title_insrt, atoms_insrt.nr, atoms_insrt.nres);
srenew(atoms_insrt.resinfo, atoms_insrt.nres);
- /* initialise van der waals arrays for inserted molecules */
- mk_vdw(&atoms_insrt, r_insrt, aps, r_distance, r_scale);
+ /* initialise distance arrays for inserted molecules */
+ exclusionDistances_insrt = makeExclusionDistances(&atoms_insrt, aps, defaultDistance, scaleFactor);
/* With -ip, take nmol_insrt from file posfn */
if (posfn != NULL)
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));
+ srenew(*exclusionDistances, (atoms->nr+atoms_insrt.nr*nmol_insrt));
trial = mol = 0;
while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
* 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))
+ if (bCheckAllPairDist && !allPairsDistOk(atoms, *x, *exclusionDistances, ePBC, box, &atoms_insrt, x_n, exclusionDistances_insrt))
{
continue;
}
- add_conf(atoms, x, NULL, r, FALSE, ePBC, box, TRUE,
- &atoms_insrt, x_n, NULL, r_insrt, FALSE, rshell, 0, oenv);
+ add_conf(atoms, x, NULL, exclusionDistances, FALSE, ePBC, box, TRUE,
+ &atoms_insrt, x_n, NULL, exclusionDistances_insrt, FALSE, rshell, 0, oenv);
if (atoms->nr == (atoms_insrt.nr+onr))
{
srenew(atoms->atomname, atoms->nr);
srenew(atoms->atom, atoms->nr);
srenew(*x, atoms->nr);
- srenew(*r, atoms->nr);
+ srenew(*exclusionDistances, 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);
+ sfree(x_n);
+ sfree(exclusionDistances_insrt);
+ done_atom(&atoms_insrt);
+
return title_insrt;
}
"By default, the insertion positions are random (with initial seed",
"specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
"molecules have been inserted in the box. Molecules are not inserted",
- "where the distance between any existing atom and any atom of ",
- "the inserted 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].[PAR]",
+ "where the distance between any existing atom and any atom of the",
+ "inserted molecule is less than the sum based on 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 the resulting radii scaled",
+ "by [TT]-scale[tt]. If radii are not found in the database, those"
+ "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
"A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
"before giving up. Increase [TT]-try[tt] if you have several small",
/* parameter data */
gmx_bool bProt, bBox;
const char *conf_prot, *confout;
- real *r;
- char *title_ins = NULL;
+ real *exclusionDistances = NULL;
+ char *title_ins = NULL;
gmx_atomprop_t aps;
/* protein configuration data */
char *title = NULL;
- t_atoms atoms;
+ t_atoms *atoms;
rvec *x = NULL;
int ePBC = -1;
matrix box;
};
#define NFILE asize(fnm)
- static int nmol_ins = 0, nmol_try = 10, seed = 1997, enum_rot;
- static real r_distance = 0.105, r_scale = 0.57;
- static rvec new_box = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
- static gmx_bool bCheckAllPairDist = FALSE;
+ static int nmol_ins = 0, nmol_try = 10, seed = 1997, enum_rot;
+ static real defaultDistance = 0.105, scaleFactor = 0.57;
+ static rvec new_box = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
+ static gmx_bool bCheckAllPairDist = FALSE;
output_env_t oenv;
const char *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
t_pargs pa[] = {
"Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
{ "-seed", FALSE, etINT, {&seed},
"Random generator seed"},
- { "-vdwd", FALSE, etREAL, {&r_distance},
+ { "-radius", FALSE, etREAL, {&defaultDistance},
"Default van der Waals distance"},
- { "-vdwscale", FALSE, etREAL, {&r_scale},
- "HIDDENScale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
+ { "-scale", FALSE, etREAL, {&scaleFactor},
+ "Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
{ "-dr", FALSE, etRVEC, {deltaR},
"Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
{ "-rot", FALSE, etENUM, {enum_rot_string},
aps = gmx_atomprop_init();
+ snew(atoms, 1);
+ init_t_atoms(atoms, 0, FALSE);
if (bProt)
{
- /*generate a solute configuration */
+ /* Generate a solute configuration */
conf_prot = opt2fn("-f", NFILE, fnm);
- title = read_conformation(conf_prot, &atoms, &x, NULL, &r, &ePBC, box,
- aps, r_distance, r_scale);
- if (atoms.nr == 0)
+ title = readConformation(conf_prot, atoms, &x, NULL,
+ &ePBC, box);
+ exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
+ if (atoms->nr == 0)
{
fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
bProt = FALSE;
}
}
- else
- {
- 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;
/* add nmol_ins molecules of atoms_ins
in random orientation at random place */
title_ins = insert_mols(insertionMoleculeFileName, nmol_ins, nmol_try, seed,
- &atoms, &x, &r, ePBC, box, aps,
- r_distance, r_scale, 0,
+ atoms, &x, &exclusionDistances, ePBC, box, aps,
+ defaultDistance, scaleFactor, 0,
oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
bCheckAllPairDist);
fprintf(stderr, "Writing generated configuration to %s\n", confout);
if (bProt)
{
- write_sto_conf(confout, title, &atoms, x, NULL, ePBC, box);
+ write_sto_conf(confout, title, atoms, x, NULL, ePBC, box);
/* print box sizes and box type to stderr */
fprintf(stderr, "%s\n", title);
}
else
{
- write_sto_conf(confout, title_ins, &atoms, x, NULL, ePBC, box);
+ write_sto_conf(confout, title_ins, atoms, x, NULL, ePBC, box);
}
/* print size of generated configuration */
fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
- atoms.nr, atoms.nres);
+ atoms->nr, atoms->nres);
gmx_atomprop_destroy(aps);
+ sfree(exclusionDistances);
+ sfree(x);
+ done_atom(atoms);
+ sfree(atoms);
return 0;
}
#include "types/atoms.h"
#include "smalloc.h"
-void mk_vdw(t_atoms *a, real rvdw[], gmx_atomprop_t aps,
- real r_distance, real r_scale)
+real *makeExclusionDistances(const t_atoms *a, gmx_atomprop_t aps,
+ real defaultDistance, real scaleFactor)
{
- int i;
+ int i;
+ real *exclusionDistances;
- /* initialise van der waals arrays of configuration */
- fprintf(stderr, "Initialising van der waals distances...\n");
+ snew(exclusionDistances, a->nr);
+ /* initialise arrays with distances usually based on van der Waals
+ radii */
+ fprintf(stderr, "Initialising inter-atomic 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])))
+ *(a->atomname[i]), &(exclusionDistances[i])))
{
- rvdw[i] = r_distance;
+ exclusionDistances[i] = defaultDistance;
}
else
{
- rvdw[i] *= r_scale;
+ exclusionDistances[i] *= scaleFactor;
}
}
+ return exclusionDistances;
}
-char *read_conformation(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
- real **r, int *ePBC, matrix box, gmx_atomprop_t aps,
- real r_distance, real r_scale)
+char *readConformation(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
+ int *ePBC, matrix box)
{
char *title;
int natoms;
{
snew(*v, natoms);
}
- snew(*r, natoms);
init_t_atoms(atoms, natoms, FALSE);
/* read residue number, residue names, atomnames, coordinates etc. */
fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
title, atoms->nr, atoms->nres);
- mk_vdw(atoms, *r, aps, r_distance, r_scale);
-
return title;
}
extern "C" {
#endif
-/*! Helper function reading VDW radii
+/*! \brief Allocate and fill an array of inter-atomic half distances
*
- * Used directly and indirectly by generate-velocities and
- * insert-molecules. */
-void mk_vdw(t_atoms *a, real rvdw[], gmx_atomprop_t aps,
- real r_distance, real r_scale);
+ * These are either scaled VDW radii taken from vdwradii.dat, or the
+ * default value. Used directly and indirectly by solvate and
+ * insert-molecules for deciding whether molecules clash. The return
+ * pointer should be freed by the caller. */
+real *makeExclusionDistances(const t_atoms *a, gmx_atomprop_t aps,
+ real defaultDistance, real scaleFactor);
-/*! Helper function to read a conformation from a file.
+/*! \brief Read a conformation from a file, allocate and fill data structures.
*
- * Used by generate-velocities and insert-molecules. */
-char *read_conformation(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
- real **r, int *ePBC, matrix box, gmx_atomprop_t aps,
- real r_distance, real r_scale);
+ * Used by solvate and insert-molecules. The returned pointers *x and
+ * *v should be freed by the caller. atoms should have its destructor
+ * called. */
+char *readConformation(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
+ int *ePBC, matrix box);
#ifdef __cplusplus
}
}
}
-static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **r,
+static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **exclusionDistances,
int ePBC, matrix box,
gmx_atomprop_t aps,
- real r_distance, real r_scale, int *atoms_added,
+ real defaultDistance, real scaleFactor, int *atoms_added,
int *residues_added, real rshell, int max_sol,
const output_env_t oenv)
{
char title_solvt[STRLEN];
t_atoms *atoms_solvt;
rvec *x_solvt, *v_solvt = NULL;
- real *r_solvt;
+ real *exclusionDistances_solvt;
int ePBC_solvt;
matrix box_solvt;
int onr, onres;
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);
+ int natoms;
+ get_stx_coordnum(filename, &natoms);
+ if (0 == natoms)
+ {
+ gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
+ }
+ snew(atoms_solvt, 1);
+ init_t_atoms(atoms_solvt, natoms, FALSE);
}
snew(x_solvt, atoms_solvt->nr);
if (v)
{
snew(v_solvt, atoms_solvt->nr);
}
- snew(r_solvt, atoms_solvt->nr);
+ snew(exclusionDistances_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);
/* apply pbc for solvent configuration for whole molecules */
rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
- /* initialise van der waals arrays for solvent configuration */
- mk_vdw(atoms_solvt, r_solvt, aps, r_distance, r_scale);
+ /* initialise distance arrays for solvent configuration */
+ exclusionDistances_solvt = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);
/* calculate the box multiplication factors n_box[0...DIM] */
nmol = 1;
{
srenew(v_solvt, atoms_solvt->nr*nmol);
}
- srenew(r_solvt, atoms_solvt->nr*nmol);
+ srenew(exclusionDistances_solvt, atoms_solvt->nr*nmol);
/* generate a new solvent configuration */
- make_new_conformation(atoms_solvt, x_solvt, v_solvt, r_solvt, box_solvt, n_box);
+ make_new_conformation(atoms_solvt, x_solvt, v_solvt, exclusionDistances_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
/* Sort the solvent mixture, not the protein... */
- sort_molecule(&atoms_solvt, x_solvt, v_solvt, r_solvt);
+ sort_molecule(&atoms_solvt, x_solvt, v_solvt, exclusionDistances_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);
+ add_conf(atoms, x, v, exclusionDistances, TRUE, ePBC, box, FALSE,
+ atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, TRUE, rshell, max_sol, oenv);
*atoms_added = atoms->nr-onr;
*residues_added = atoms->nres-onres;
sfree(x_solvt);
- sfree(r_solvt);
+ sfree(exclusionDistances_solvt);
+ done_atom(atoms_solvt);
+ sfree(atoms_solvt);
fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
*atoms_added, *residues_added);
"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]",
+ "the solvent molecule is less than the sum of the scaled van der Waals",
+ "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
+ "Waals radii is read by the program, and the resulting radii scaled",
+ "by [TT]-scale[tt]. If radii are not found in the database, those"
+ "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
"The default solvent is Simple Point Charge water (SPC), with coordinates ",
"from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
/* parameter data */
gmx_bool bProt, bBox;
const char *conf_prot, *confout;
- real *r;
+ real *exclusionDistances = NULL;
gmx_atomprop_t aps;
/* protein configuration data */
char *title = NULL;
- t_atoms atoms;
- rvec *x, *v = NULL;
+ t_atoms *atoms;
+ rvec *x = NULL, *v = NULL;
int ePBC = -1;
matrix box;
};
#define NFILE asize(fnm)
- static real r_distance = 0.105, r_shell = 0, r_scale = 0.57;
- static rvec new_box = {0.0, 0.0, 0.0};
- static gmx_bool bReadV = FALSE;
- static int max_sol = 0;
+ static real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
+ static rvec new_box = {0.0, 0.0, 0.0};
+ static gmx_bool bReadV = FALSE;
+ static int max_sol = 0;
output_env_t oenv;
t_pargs pa[] = {
{ "-box", FALSE, etRVEC, {new_box},
"Box size (in nm)" },
- { "-vdwd", FALSE, etREAL, {&r_distance},
+ { "-radius", FALSE, etREAL, {&defaultDistance},
"Default van der Waals distance"},
- { "-vdwscale", FALSE, etREAL, {&r_scale},
- "HIDDENScale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
+ { "-scale", FALSE, etREAL, {&scaleFactor},
+ "Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
{ "-shell", FALSE, etREAL, {&r_shell},
"Thickness of optional water layer around solute" },
{ "-maxsol", FALSE, etINT, {&max_sol},
aps = gmx_atomprop_init();
+ snew(atoms, 1);
+ init_t_atoms(atoms, 0, FALSE);
if (bProt)
{
/* Generate a solute configuration */
conf_prot = opt2fn("-cp", NFILE, fnm);
- title = read_conformation(conf_prot, &atoms, &x, bReadV ? &v : NULL, &r, &ePBC, box,
- aps, r_distance, r_scale);
+ title = readConformation(conf_prot, atoms, &x,
+ bReadV ? &v : NULL, &ePBC, box);
+ exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
+
if (bReadV && !v)
{
fprintf(stderr, "Note: no velocities found\n");
}
- if (atoms.nr == 0)
+ if (atoms->nr == 0)
{
fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
bProt = FALSE;
}
}
- else
- {
- 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;
"or give explicit -box command line option");
}
- add_solv(solventFileName, &atoms, &x, v ? &v : NULL, &r, ePBC, box,
- aps, r_distance, r_scale, &atoms_added, &residues_added, r_shell, max_sol,
+ add_solv(solventFileName, atoms, &x, v ? &v : NULL, &exclusionDistances, ePBC, box,
+ aps, defaultDistance, scaleFactor, &atoms_added, &residues_added, r_shell, max_sol,
oenv);
/* write new configuration 1 to file confout */
fprintf(stderr, "Writing generated configuration to %s\n", confout);
if (bProt)
{
- write_sto_conf(confout, title, &atoms, x, v, ePBC, box);
+ 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, "Generated by gmx solvate", &atoms, x, v, ePBC, box);
+ write_sto_conf(confout, "Generated by gmx solvate", atoms, x, v, ePBC, box);
}
/* 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);
+ atoms->nr, atoms->nres);
+ update_top(atoms, box, NFILE, fnm, aps);
gmx_atomprop_destroy(aps);
+ sfree(exclusionDistances);
+ sfree(x);
+ sfree(v);
+ done_atom(atoms);
+ sfree(atoms);
return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff