2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2007, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
72 void add_rec(t_conect c[],atom_id i,atom_id j,real d2)
74 if (c[i].aa == NO_ATID) {
78 else if (c[i].ab == NO_ATID) {
82 else if (d2 < c[i].d2a) {
86 else if (d2 < c[i].d2b) {
90 /* Swap them if necessary: a must be larger than b */
91 if (c[i].d2a < c[i].d2b) {
101 void do_conect(const char *fn,int n,rvec x[])
109 fprintf(stderr,"Building CONECT records\n");
112 c[i].aa = c[i].ab = NO_ATID;
114 for(i=0; (i<n); i++) {
115 for(j=i+1; (j<n); j++) {
116 rvec_sub(x[i],x[j],dx);
123 for(i=0; (i<n); i++) {
124 if ((c[i].aa == NO_ATID) || (c[i].ab == NO_ATID))
125 fprintf(stderr,"Warning dot %d has no conections\n",i+1);
126 fprintf(fp,"CONECT%5d%5d%5d\n",i+1,c[i].aa+1,c[i].ab+1);
132 void connelly_plot(const char *fn,int ndots,real dots[],rvec x[],t_atoms *atoms,
133 t_symtab *symtab,int ePBC,matrix box,gmx_bool bSave)
135 static const char *atomnm="DOT";
136 static const char *resnm ="DOT";
137 static const char *title ="Connely Dot Surface Generated by g_sas";
146 srenew(atoms->atom,atoms->nr+ndots);
147 srenew(atoms->atomname,atoms->nr+ndots);
148 srenew(atoms->resinfo,r0+1);
149 atoms->atom[i0].resind = r0;
150 t_atoms_set_resinfo(atoms,i0,symtab,resnm,r0+1,' ',0,' ');
151 srenew(atoms->pdbinfo,atoms->nr+ndots);
152 snew(xnew,atoms->nr+ndots);
153 for(i=0; (i<atoms->nr); i++)
154 copy_rvec(x[i],xnew[i]);
155 for(i=k=0; (i<ndots); i++) {
157 atoms->atomname[ii0] = put_symtab(symtab,atomnm);
158 atoms->pdbinfo[ii0].type = epdbATOM;
159 atoms->pdbinfo[ii0].atomnr= ii0;
160 atoms->atom[ii0].resind = r0;
161 xnew[ii0][XX] = dots[k++];
162 xnew[ii0][YY] = dots[k++];
163 xnew[ii0][ZZ] = dots[k++];
164 atoms->pdbinfo[ii0].bfac = 0.0;
165 atoms->pdbinfo[ii0].occup = 0.0;
167 atoms->nr = i0+ndots;
169 write_sto_conf(fn,title,atoms,xnew,NULL,ePBC,box);
174 init_t_atoms(&aaa,ndots,TRUE);
175 aaa.atom[0].resind = 0;
176 t_atoms_set_resinfo(&aaa,0,symtab,resnm,1,' ',0,' ');
178 for(i=k=0; (i<ndots); i++) {
180 aaa.atomname[ii0] = put_symtab(symtab,atomnm);
181 aaa.pdbinfo[ii0].type = epdbATOM;
182 aaa.pdbinfo[ii0].atomnr= ii0;
183 aaa.atom[ii0].resind = 0;
184 xnew[ii0][XX] = dots[k++];
185 xnew[ii0][YY] = dots[k++];
186 xnew[ii0][ZZ] = dots[k++];
187 aaa.pdbinfo[ii0].bfac = 0.0;
188 aaa.pdbinfo[ii0].occup = 0.0;
191 write_sto_conf(fn,title,&aaa,xnew,NULL,ePBC,box);
192 do_conect(fn,ndots,xnew);
193 free_t_atoms(&aaa,FALSE);
198 real calc_radius(char *atom)
224 void sas_plot(int nfile,t_filenm fnm[],real solsize,int ndots,
225 real qcut,gmx_bool bSave,real minarea,gmx_bool bPBC,
226 real dgs_default,gmx_bool bFindex, const output_env_t oenv)
228 FILE *fp,*fp2,*fp3=NULL,*vp;
229 const char *flegend[] = { "Hydrophobic", "Hydrophilic",
230 "Total", "D Gsolv" };
231 const char *vlegend[] = { "Volume (nm\\S3\\N)", "Density (g/l)" };
232 const char *or_and_oa_legend[] = { "Average (nm\\S2\\N)", "Standard deviation (nm\\S2\\N)" };
235 gmx_atomprop_t aps=NULL;
236 gmx_rmpbc_t gpbc=NULL;
239 int i,j,ii,nfr,natoms,flag,nsurfacedots,res;
247 gmx_bool *bOut,*bPhobic;
249 gmx_bool bResAt,bITP,bDGsol;
250 real *radius,*dgs_factor=NULL,*area=NULL,*surfacedots=NULL;
251 real at_area,*atom_area=NULL,*atom_area2=NULL;
252 real *res_a=NULL,*res_area=NULL,*res_area2=NULL;
253 real totarea,totvolume,totmass=0,density,harea,tarea,fluc2;
254 atom_id **index,*findex;
255 int *nx,nphobic,npcheck,retval;
256 char **grpname,*fgrpname;
259 bITP = opt2bSet("-i",nfile,fnm);
260 bResAt = opt2bSet("-or",nfile,fnm) || opt2bSet("-oa",nfile,fnm) || bITP;
262 bTop = read_tps_conf(ftp2fn(efTPS,nfile,fnm),title,&top,&ePBC,
263 &xtop,NULL,topbox,FALSE);
264 atoms = &(top.atoms);
267 fprintf(stderr,"No tpr file, will not compute Delta G of solvation\n");
270 bDGsol = strcmp(*(atoms->atomtype[0]),"?") != 0;
272 fprintf(stderr,"Warning: your tpr file is too old, will not compute "
273 "Delta G of solvation\n");
275 printf("In case you use free energy of solvation predictions:\n");
276 please_cite(stdout,"Eisenberg86a");
280 aps = gmx_atomprop_init();
282 if ((natoms=read_first_x(oenv,&status,ftp2fn(efTRX,nfile,fnm),
284 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
286 if ((ePBC != epbcXYZ) || (TRICLINIC(box))) {
287 fprintf(stderr,"\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
288 "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
289 "will certainly crash the analysis.\n\n");
294 fprintf(stderr,"Select a group for calculation of surface and a group for output:\n");
295 get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),2,nx,index,grpname);
298 fprintf(stderr,"Select a group of hydrophobic atoms:\n");
299 get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),1,&nphobic,&findex,&fgrpname);
302 for(i=0; i<nx[1]; i++)
303 bOut[index[1][i]] = TRUE;
305 /* Now compute atomic readii including solvent probe size */
309 snew(atom_area,nx[0]);
310 snew(atom_area2,nx[0]);
311 snew(res_a,atoms->nres);
312 snew(res_area,atoms->nres);
313 snew(res_area2,atoms->nres);
316 snew(dgs_factor,nx[0]);
318 /* Get a Van der Waals radius for each atom */
320 for(i=0; (i<natoms); i++) {
321 if (!gmx_atomprop_query(aps,epropVDW,
322 *(atoms->resinfo[atoms->atom[i].resind].name),
323 *(atoms->atomname[i]),&radius[i]))
325 /* radius[i] = calc_radius(*(top->atoms.atomname[i])); */
326 radius[i] += solsize;
329 fprintf(stderr,"WARNING: could not find a Van der Waals radius for %d atoms\n",ndefault);
330 /* Determine which atom is counted as hydrophobic */
333 for(i=0; (i<nx[0]); i++) {
335 for(j=0; (j<nphobic); j++) {
336 if (findex[j] == ii) {
343 if (npcheck != nphobic)
344 gmx_fatal(FARGS,"Consistency check failed: not all %d atoms in the hydrophobic index\n"
345 "found in the normal index selection (%d atoms)",nphobic,npcheck);
350 for(i=0; (i<nx[0]); i++) {
353 bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
354 if (bPhobic[i] && bOut[ii])
358 if (!gmx_atomprop_query(aps,epropDGsol,
359 *(atoms->resinfo[atoms->atom[ii].resind].name),
360 *(atoms->atomtype[ii]),&(dgs_factor[i])))
361 dgs_factor[i] = dgs_default;
363 fprintf(debug,"Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
364 ii+1,*(atoms->resinfo[atoms->atom[ii].resind].name),
365 *(atoms->atomname[ii]),
366 atoms->atom[ii].q,radius[ii]-solsize,dgs_factor[i],
369 fprintf(stderr,"%d out of %d atoms were classified as hydrophobic\n",
372 fp=xvgropen(opt2fn("-o",nfile,fnm),"Solvent Accessible Surface","Time (ps)",
373 "Area (nm\\S2\\N)",oenv);
374 xvgr_legend(fp,asize(flegend) - (bDGsol ? 0 : 1),flegend,oenv);
375 vfile = opt2fn_null("-tv",nfile,fnm);
378 gmx_fatal(FARGS,"Need a tpr file for option -tv");
380 vp=xvgropen(vfile,"Volume and Density","Time (ps)","",oenv);
381 xvgr_legend(vp,asize(vlegend),vlegend,oenv);
384 for(i=0; (i<nx[0]); i++) {
388 if (!query_atomprop(atomprop,epropMass,
389 *(top->atoms.resname[top->atoms.atom[ii].resnr]),
390 *(top->atoms.atomname[ii]),&mm))
394 totmass += atoms->atom[ii].m;
397 fprintf(stderr,"WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n",ndefault);
402 gmx_atomprop_destroy(aps);
405 gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
410 gmx_rmpbc(gpbc,natoms,box,x);
412 bConnelly = (nfr==0 && opt2bSet("-q",nfile,fnm));
415 gmx_fatal(FARGS,"Need a tpr file for Connelly plot");
416 flag = FLAG_ATOM_AREA | FLAG_DOTS;
418 flag = FLAG_ATOM_AREA;
421 flag = flag | FLAG_VOLUME;
425 write_sto_conf("check.pdb","pbc check",atoms,x,NULL,ePBC,box);
427 retval = nsc_dclm_pbc(x,radius,nx[0],ndots,flag,&totarea,
428 &area,&totvolume,&surfacedots,&nsurfacedots,
429 index[0],ePBC,bPBC ? box : NULL);
431 gmx_fatal(FARGS,"Something wrong in nsc_dclm_pbc");
434 connelly_plot(ftp2fn(efPDB,nfile,fnm),
435 nsurfacedots,surfacedots,x,atoms,
436 &(top.symtab),ePBC,box,bSave);
441 for(i=0; i<atoms->nres; i++)
443 for(i=0; (i<nx[0]); i++) {
448 atom_area[i] += at_area;
449 atom_area2[i] += sqr(at_area);
450 res_a[atoms->atom[ii].resind] += at_area;
454 dgsolv += at_area*dgs_factor[i];
460 for(i=0; i<atoms->nres; i++) {
461 res_area[i] += res_a[i];
462 res_area2[i] += sqr(res_a[i]);
464 fprintf(fp,"%10g %10g %10g %10g",t,harea,tarea-harea,tarea);
466 fprintf(fp," %10g\n",dgsolv);
472 density = totmass*AMU/(totvolume*NANO*NANO*NANO);
473 fprintf(vp,"%12.5e %12.5e %12.5e\n",t,totvolume,density);
484 } while (read_next_x(oenv,status,&t,natoms,x,box));
487 gmx_rmpbc_done(gpbc);
489 fprintf(stderr,"\n");
495 /* if necessary, print areas per atom to file too: */
497 for(i=0; i<atoms->nres; i++) {
501 for(i=0; i<nx[0]; i++) {
503 atom_area2[i] /= nfr;
505 fprintf(stderr,"Printing out areas per atom\n");
506 fp = xvgropen(opt2fn("-or",nfile,fnm),"Area per residue over the trajectory","Residue",
507 "Area (nm\\S2\\N)",oenv);
508 xvgr_legend(fp, asize(or_and_oa_legend),or_and_oa_legend,oenv);
509 fp2 = xvgropen(opt2fn("-oa",nfile,fnm),"Area per atom over the trajectory","Atom #",
510 "Area (nm\\S2\\N)",oenv);
511 xvgr_legend(fp2, asize(or_and_oa_legend),or_and_oa_legend,oenv);
513 fp3 = ftp2FILE(efITP,nfile,fnm,"w");
514 fprintf(fp3,"[ position_restraints ]\n"
518 "; Atom Type fx fy fz\n");
520 for(i=0; i<nx[0]; i++) {
522 res = atoms->atom[ii].resind;
523 if (i==nx[0]-1 || res!=atoms->atom[index[0][i+1]].resind) {
524 fluc2 = res_area2[res]-sqr(res_area[res]);
527 fprintf(fp,"%10d %10g %10g\n",
528 atoms->resinfo[res].nr,res_area[res],sqrt(fluc2));
530 fluc2 = atom_area2[i]-sqr(atom_area[i]);
533 fprintf(fp2,"%d %g %g\n",index[0][i]+1,atom_area[i],sqrt(fluc2));
534 if (bITP && (atom_area[i] > minarea))
535 fprintf(fp3,"%5d 1 FCX FCX FCZ\n",ii+1);
542 /* Be a good citizen, keep our memory free! */
568 int gmx_sas(int argc,char *argv[])
570 const char *desc[] = {
571 "[TT]g_sas[tt] computes hydrophobic, hydrophilic and total solvent",
572 "accessible surface area. See Eisenhaber F, Lijnzaad P, Argos P,",
573 "Sander C, & Scharf M (1995) J. Comput. Chem. 16, 273-284.",
574 "As a side effect, the Connolly surface can be generated as well in",
575 "a [TT].pdb[tt] file where the nodes are represented as atoms and the",
576 "vertice connecting the nearest nodes as CONECT records.",
577 "The program will ask for a group for the surface calculation",
578 "and a group for the output. The calculation group should always",
579 "consists of all the non-solvent atoms in the system.",
580 "The output group can be the whole or part of the calculation group.",
581 "The average and standard deviation of the area over the trajectory can be plotted",
582 "per residue and atom as well (options [TT]-or[tt] and [TT]-oa[tt]).",
583 "In combination with the latter option an [TT].itp[tt] file can be",
584 "generated (option [TT]-i[tt])",
585 "which can be used to restrain surface atoms.[PAR]",
587 "By default, periodic boundary conditions are taken into account,",
588 "this can be turned off using the [TT]-nopbc[tt] option.[PAR]",
590 "With the [TT]-tv[tt] option the total volume and density of the",
591 "molecule can be computed.",
592 "Please consider whether the normal probe radius is appropriate",
593 "in this case or whether you would rather use e.g. 0. It is good",
594 "to keep in mind that the results for volume and density are very",
595 "approximate. For example, in ice Ih, one can easily fit water molecules in the",
596 "pores which would yield a volume that is too low, and surface area and density",
597 "that are both too high."
601 static real solsize = 0.14;
602 static int ndots = 24;
603 static real qcut = 0.2;
604 static real minarea = 0.5, dgs_default=0;
605 static gmx_bool bSave = TRUE,bPBC=TRUE,bFindex=FALSE;
607 { "-probe", FALSE, etREAL, {&solsize},
608 "Radius of the solvent probe (nm)" },
609 { "-ndots", FALSE, etINT, {&ndots},
610 "Number of dots per sphere, more dots means more accuracy" },
611 { "-qmax", FALSE, etREAL, {&qcut},
612 "The maximum charge (e, absolute value) of a hydrophobic atom" },
613 { "-f_index", FALSE, etBOOL, {&bFindex},
614 "Determine from a group in the index file what are the hydrophobic atoms rather than from the charge" },
615 { "-minarea", FALSE, etREAL, {&minarea},
616 "The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)" },
617 { "-pbc", FALSE, etBOOL, {&bPBC},
618 "Take periodicity into account" },
619 { "-prot", FALSE, etBOOL, {&bSave},
620 "Output the protein to the Connelly [TT].pdb[tt] file too" },
621 { "-dgs", FALSE, etREAL, {&dgs_default},
622 "Default value for solvation free energy per area (kJ/mol/nm^2)" }
625 { efTRX, "-f", NULL, ffREAD },
626 { efTPS, "-s", NULL, ffREAD },
627 { efXVG, "-o", "area", ffWRITE },
628 { efXVG, "-or", "resarea", ffOPTWR },
629 { efXVG, "-oa", "atomarea", ffOPTWR },
630 { efXVG, "-tv", "volume", ffOPTWR },
631 { efPDB, "-q", "connelly", ffOPTWR },
632 { efNDX, "-n", "index", ffOPTRD },
633 { efITP, "-i", "surfat", ffOPTWR }
635 #define NFILE asize(fnm)
637 CopyRight(stderr,argv[0]);
638 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
639 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
642 fprintf(stderr,"Probe size too small, setting it to %g\n",solsize);
646 fprintf(stderr,"Ndots too small, setting it to %d\n",ndots);
649 please_cite(stderr,"Eisenhaber95");
651 sas_plot(NFILE,fnm,solsize,ndots,qcut,bSave,minarea,bPBC,dgs_default,bFindex,
654 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
655 do_view(oenv,opt2fn_null("-or",NFILE,fnm),"-nxy");
656 do_view(oenv,opt2fn_null("-oa",NFILE,fnm),"-nxy");