int res0;
} t_moltypes;
-static void sort_molecule(t_atoms **atoms_solvt, std::vector<RVec> *x,
+static void sort_molecule(t_atoms **atoms_solvt,
+ t_atoms **newatoms,
+ std::vector<RVec> *x,
std::vector<RVec> *v)
{
int atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
t_moltypes *moltypes;
- t_atoms *atoms, *newatoms;
+ t_atoms *atoms;
fprintf(stderr, "Sorting configuration\n");
}
/* now put them there: */
- snew(newatoms, 1);
- init_t_atoms(newatoms, atoms->nr, FALSE);
- newatoms->nres = atoms->nres;
- snew(newatoms->resinfo, atoms->nres);
+ snew(*newatoms, 1);
+ init_t_atoms(*newatoms, atoms->nr, FALSE);
+ (*newatoms)->nres = atoms->nres;
+ srenew((*newatoms)->resinfo, atoms->nres);
std::vector<RVec> newx(x->size());
std::vector<RVec> newv(v->size());
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;
+ (*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;
+ (*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->empty())
{
/* put them back into the original arrays and throw away temporary arrays */
done_atom(atoms);
- *atoms_solvt = newatoms;
+ *atoms_solvt = (*newatoms);
std::swap(*x, newx);
std::swap(*v, newv);
}
}
/* Sort the solvent mixture, not the protein... */
- sort_molecule(&atoms_solvt, &x_solvt, &v_solvt);
+ t_atoms *newatoms = nullptr;
+ sort_molecule(&atoms_solvt, &newatoms, &x_solvt, &v_solvt);
// Merge the two configurations.
x->insert(x->end(), x_solvt.begin(), x_solvt.end());
done_top(top_solvt);
sfree(top_solvt);
+ if (newatoms)
+ {
+ done_atom(newatoms);
+ sfree(newatoms);
+ }
}
-static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
+static void update_top(t_atoms *atoms, int firstSolventResidueIndex, matrix box, int NFILE, t_filenm fnm[],
gmx_atomprop_t aps)
{
- FILE *fpin, *fpout;
- char buf[STRLEN], buf2[STRLEN], *temp;
- const char *topinout;
- int line;
- 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++;
- }
- }
+ FILE *fpin, *fpout;
+ char buf[STRLEN], buf2[STRLEN], *temp;
+ const char *topinout;
+ int line;
+ bool bSystem;
+ int i;
+ double mtot;
+ real vol, mm;
+
+ int nsol = atoms->nres - firstSolventResidueIndex;
+
mtot = 0;
for (i = 0; (i < atoms->nr); i++)
{
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);
+ fprintf(stderr, "Number of solvent molecules: %5d \n\n", nsol);
/* open topology file and append sol molecules */
topinout = ftp2fn(efTOP, NFILE, fnm);
fpin = gmx_ffopen(topinout, "r");
fpout = gmx_fopen_temporary(temporary_filename);
line = 0;
- bSystem = bMolecules = false;
+ bSystem = false;
while (fgets(buf, STRLEN, fpin))
{
- bSkip = false;
line++;
strcpy(buf2, buf);
if ((temp = strchr(buf2, '\n')) != nullptr)
ltrim(buf2);
rtrim(buf2);
bSystem = (gmx_strcasecmp(buf2, "system") == 0);
- bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
}
}
else if (bSystem && nsol && (buf[0] != ';') )
bSystem = false;
}
}
- else if (bMolecules)
+ fprintf(fpout, "%s", buf);
+ }
+ gmx_ffclose(fpin);
+
+ // Add new solvent molecules to the topology
+ if (nsol > 0)
+ {
+ std::string currRes = *atoms->resinfo[firstSolventResidueIndex].name;
+ int resCount = 0;
+
+ // Iterate through solvent molecules and increment a count until new resname found
+ for (int i = firstSolventResidueIndex; i < atoms->nres; i++)
{
- /* check if this is a line with solvent molecules */
- sscanf(buf, "%4095s", buf2);
- if (strcmp(buf2, "SOL") == 0)
+ if ((currRes == *atoms->resinfo[i].name))
{
- sscanf(buf, "%*4095s %20d", &i);
- nsol -= i;
- if (nsol < 0)
- {
- bSkip = true;
- nsol += i;
- }
+ resCount += 1;
}
- }
- if (bSkip)
- {
- if ((temp = strchr(buf, '\n')) != nullptr)
+ else
{
- temp[0] = '\0';
+ // Change topology and restart count
+ fprintf(stdout, "Adding line for %d solvent molecules with resname (%s) to "
+ "topology file (%s)\n", resCount, currRes.c_str(), topinout);
+ fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
+ currRes = *atoms->resinfo[i].name;
+ resCount = 1;
}
- fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
- line, buf, topinout);
- }
- else
- {
- fprintf(fpout, "%s", buf);
}
- }
- gmx_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);
+ // One more print needed for last residue type
+ fprintf(stdout, "Adding line for %d solvent molecules with resname (%s) to "
+ "topology file (%s)\n", resCount, currRes.c_str(), topinout);
+ fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
}
gmx_ffclose(fpout);
make_backup(topinout);
};
#define NFILE asize(fnm)
- real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
- rvec new_box = {0.0, 0.0, 0.0};
- gmx_bool bReadV = FALSE;
- int max_sol = 0;
+ real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
+ rvec new_box = {0.0, 0.0, 0.0};
+ gmx_bool bReadV = FALSE;
+ int max_sol = 0;
+ int firstSolventResidueIndex = 0;
gmx_output_env_t *oenv;
t_pargs pa[] = {
{ "-box", FALSE, etRVEC, {new_box},
fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
bProt = FALSE;
}
+ else
+ {
+ firstSolventResidueIndex = top->atoms.nres;
+ }
}
if (bBox)
{
/* print size of generated configuration */
fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
top->atoms.nr, top->atoms.nres);
- update_top(&top->atoms, box, NFILE, fnm, aps);
+ update_top(&top->atoms, firstSolventResidueIndex, box, NFILE, fnm, aps);
gmx_atomprop_destroy(aps);
done_top(top);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff