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-2004, 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.
55 #include "gmx_fatal.h"
66 void calc_rm_cm(int isize, atom_id index[], t_atoms *atoms, rvec x[], rvec xcm)
71 /* calculate and remove center of mass of reference structure */
74 for(i=0; i<isize; i++) {
75 m = atoms->atom[index[i]].m;
77 xcm[d] += m*x[index[i]][d];
81 for(i=0; i<atoms->nr; i++)
85 int build_res_index(int isize, atom_id index[], t_atom atom[], int rindex[])
90 rindex[r] = atom[index[0]].resind;
92 for(i=1; i<isize; i++)
93 if (atom[index[i]].resind != rindex[r-1]) {
94 rindex[r] = atom[index[i]].resind;
101 int find_res_end(int i, int isize, atom_id index[], t_atoms *atoms)
105 rnr = atoms->atom[index[i]].resind;
106 while(i<isize && atoms->atom[index[i]].resind==rnr)
111 int debug_strcmp(char s1[], char s2[])
113 if (debug) fprintf(debug," %s-%s",s1,s2);
114 return strcmp(s1, s2);
117 int find_next_match_atoms_in_res(int *i1,int isize1,atom_id index1[],
118 int m1,char **atnms1[],
119 int *i2,int isize2,atom_id index2[],
120 int m2,char **atnms2[])
122 int dx, dy, dmax, cmp;
127 dmax = max(m1-*i1, m2-*i2);
128 for(dx=0; dx<dmax && cmp!=0; dx++) {
129 for(dy=dx; dy<dmax && cmp!=0; dy++) {
131 if (debug) fprintf(debug, ".");
133 if ( *i1+dx<m1 && *i2+dy<m2 ) {
135 cmp = debug_strcmp(*atnms1[index1[*i1+dx]],*atnms2[index2[*i2+dy]]);
136 if (debug) fprintf(debug,"(%d %d)", *i1+dx, *i2+dy);
138 if ( cmp!=0 && *i1+dy<m1 && *i2+dx<m2 ) {
140 cmp = debug_strcmp(*atnms1[index1[*i1+dy]],*atnms2[index2[*i2+dx]]);
141 if (debug) fprintf(debug,"(%d %d)", *i1+dy, *i2+dx);
146 /* apparently, dx and dy are incremented one more time
147 as the loops terminate: we correct this here */
151 if (debug) fprintf(debug, "{%d %d}", *i1+bFW?dx:dy, *i2+bFW?dy:dx );
164 static int find_next_match_res(int *rnr1, int isize1,
165 int index1[], t_resinfo *resinfo1,
166 int *rnr2, int isize2,
167 int index2[], t_resinfo *resinfo2)
169 int dx, dy, dmax, cmp, rr1, rr2;
170 gmx_bool bFW=FALSE,bFF=FALSE;
174 while(rr1<isize1 && *rnr1 != index1[rr1])
177 while(rr2<isize2 && *rnr2 != index2[rr2])
181 dmax = max(isize1-rr1, isize2-rr2);
182 if (debug) fprintf(debug," R:%d-%d:%d-%d:%d ",
183 rr1, isize1, rr2, isize2, dmax);
184 for(dx=0; dx<dmax && cmp!=0; dx++)
185 for(dy=0; dy<=dx && cmp!=0; dy++)
188 if ( rr1+dx<isize1 && rr2+dy<isize2 ) {
190 cmp=debug_strcmp(*resinfo1[index1[rr1+dx]].name,
191 *resinfo2[index2[rr2+dy]].name);
192 if (debug) fprintf(debug,"(%d %d)", rr1+dx, rr2+dy);
194 if ( cmp!=0 && rr1+dy<isize1 && rr2+dx<isize2 ) {
196 cmp=debug_strcmp(*resinfo1[index1[rr1+dy]].name,
197 *resinfo2[index2[rr2+dx]].name);
198 if (debug) fprintf(debug,"(%d %d)", rr1+dy, rr2+dx);
200 if ( dx!=0 && cmp!=0 && rr1+dx<isize1 && rr2+dx<isize2 ) {
202 cmp=debug_strcmp(*resinfo1[index1[rr1+dx]].name,
203 *resinfo2[index2[rr2+dx]].name);
204 if (debug) fprintf(debug,"(%d %d)", rr1+dx, rr2+dx);
208 /* apparently, dx and dy are incremented one more time
209 as the loops terminate: we correct this here */
212 /* if we skipped equal on both sides, only skip one residue
213 to allow for single mutations of residues... */
215 if (debug) fprintf(debug, "%d.%d.%dX%sX%s",dx,rr1,rr2,
216 *resinfo1[index1[rr1+1]].name,
217 *resinfo2[index2[rr2+1]].name);
221 if (debug) fprintf(debug, "!");
240 int find_first_atom_in_res(int rnr, int isize, atom_id index[], t_atom atom[])
245 while(i<isize && atom[index[i]].resind!=rnr)
248 if (atom[index[i]].resind==rnr)
254 void find_matching_names(int *isize1, atom_id index1[], t_atoms *atoms1,
255 int *isize2, atom_id index2[], t_atoms *atoms2)
257 int i, i1, i2, ii1, ii2, m1, m2;
259 int r,rnr1, rnr2, prnr1, prnr2;
261 int *rindex1, *rindex2;
262 char *resnm1, *resnm2, *atnm1, *atnm2;
263 char ***atnms1,***atnms2;
264 t_resinfo *resinfo1,*resinfo2;
266 /* set some handy shortcuts */
267 resinfo1 = atoms1->resinfo;
268 atnms1 = atoms1->atomname;
269 resinfo2 = atoms2->resinfo;
270 atnms2 = atoms2->atomname;
272 /* build indexes of selected residues */
273 snew(rindex1,atoms1->nres);
274 rsize1=build_res_index(*isize1, index1, atoms1->atom, rindex1);
275 snew(rindex2,atoms2->nres);
276 rsize2=build_res_index(*isize2, index2, atoms2->atom, rindex2);
282 if (debug) fprintf(debug, "Find matching names: %d, %d\n", *isize1, *isize2);
283 while ( atcmp==0 && i1<*isize1 && i2<*isize2 ) {
285 rnr1 = atoms1->atom[index1[i1]].resind;
286 rnr2 = atoms2->atom[index2[i2]].resind;
287 if (rnr1!=prnr1 || rnr2!=prnr2) {
288 if (debug) fprintf(debug, "R: %s%d %s%d\n",
289 *resinfo1[rnr1].name,rnr1,*resinfo2[rnr2].name,rnr2);
290 rescmp=strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
292 if (debug) fprintf(debug, "comparing %d %d", i1, i2);
293 atcmp=debug_strcmp(*atnms1[index1[i1]], *atnms2[index2[i2]]);
296 if (atcmp!=0) { /* no match -> find match within residues */
297 m1 = find_res_end(i1, *isize1, index1, atoms1);
298 m2 = find_res_end(i2, *isize2, index2, atoms2);
299 if (debug) fprintf(debug, " [%d<%d %d<%d]", i1,m1, i2,m2);
300 atcmp = find_next_match_atoms_in_res(&i1,*isize1,index1,m1,atnms1,
301 &i2,*isize2,index2,m2,atnms2);
302 if (debug) fprintf(debug, " -> %d %d %s-%s", i1, i2,
303 *atnms1[index1[i1]],*atnms2[index2[i2]]);
306 if (atcmp!=0) { /* still no match -> next residue(s) */
309 rescmp = find_next_match_res(&rnr1,rsize1,rindex1,resinfo1,
310 &rnr2,rsize2,rindex2,resinfo2);
312 i1 = find_first_atom_in_res(rnr1, *isize1, index1, atoms1->atom);
314 i2 = find_first_atom_in_res(rnr2, *isize2, index2, atoms2->atom);
315 if (debug) fprintf(debug, " -> %s%d-%s%d %s%d-%s%d",
316 *resinfo1[rnr1].name,rnr1,
317 *atnms1[index1[i1]],index1[i1],
318 *resinfo2[rnr2].name,rnr2,
319 *atnms2[index2[i2]],index2[i2]);
320 m1 = find_res_end(i1, *isize1, index1, atoms1);
321 m2 = find_res_end(i2, *isize2, index2, atoms2);
322 if (debug) fprintf(debug, " [%d<%d %d<%d]", i1,m1, i2,m2);
323 atcmp = find_next_match_atoms_in_res(&i1,*isize1,index1,m1,atnms1,
324 &i2,*isize2,index2,m2,atnms2);
325 if (debug) fprintf(debug, " -> %d %d %s-%s", i1, i2,
326 *atnms1[index1[i1]],*atnms2[index2[i2]]);
329 fprintf(debug, "(%d %d): %d %d\n", ii1, ii2, atcmp, rescmp);
330 if (atcmp==0) { /* if match -> store indices */
331 index1[ii1++]=index1[i1];
332 index2[ii2++]=index2[i2];
343 gmx_fatal(FARGS, "DEATH HORROR: non-equal number of matching atoms!\n");
344 if ( ii1==i1 && ii2==i2 )
345 printf("All atoms in index groups 1 and 2 match\n");
347 if ( i1==i2 && ii1==ii2 )
348 printf("Both index groups modified from %d to %d atoms\n", i1, ii1);
351 printf("Index group 1 modified from %d to %d atoms\n",i1, ii1);
353 printf("Index group 2 modified from %d to %d atoms\n",i2, ii2);
361 int gmx_confrms(int argc,char *argv[])
363 const char *desc[] = {
364 "[TT]g_confrms[tt] computes the root mean square deviation (RMSD) of two",
365 "structures after least-squares fitting the second structure on the first one.",
366 "The two structures do NOT need to have the same number of atoms,",
367 "only the two index groups used for the fit need to be identical.",
368 "With [TT]-name[tt] only matching atom names from the selected groups",
369 "will be used for the fit and RMSD calculation. This can be useful ",
370 "when comparing mutants of a protein.",
372 "The superimposed structures are written to file. In a [TT].pdb[tt] file",
373 "the two structures will be written as separate models",
374 "(use [TT]rasmol -nmrpdb[tt]). Also in a [TT].pdb[tt] file, B-factors",
375 "calculated from the atomic MSD values can be written with [TT]-bfac[tt].",
377 static gmx_bool bOne=FALSE,bRmpbc=FALSE,bMW=TRUE,bName=FALSE,
378 bBfac=FALSE,bFit=TRUE,bLabel=FALSE;
381 { "-one", FALSE, etBOOL, {&bOne}, "Only write the fitted structure to file" },
382 { "-mw", FALSE, etBOOL, {&bMW}, "Mass-weighted fitting and RMSD" },
383 { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
384 { "-fit", FALSE, etBOOL, {&bFit},
385 "Do least squares superposition of the target structure to the reference" },
386 { "-name", FALSE,etBOOL,{&bName},
387 "Only compare matching atom names" },
388 { "-label",FALSE,etBOOL,{&bLabel},
389 "Added chain labels A for first and B for second structure"},
390 { "-bfac", FALSE,etBOOL,{&bBfac},
391 "Output B-factors from atomic MSD values" }
394 { efTPS, "-f1", "conf1.gro", ffREAD },
395 { efSTX, "-f2", "conf2", ffREAD },
396 { efSTO, "-o", "fit.pdb", ffWRITE },
397 { efNDX, "-n1" , "fit1.ndx", ffOPTRD },
398 { efNDX, "-n2", "fit2.ndx", ffOPTRD },
399 { efNDX, "-no", "match.ndx", ffOPTWR }
401 #define NFILE asize(fnm)
403 /* the two structure files */
404 const char *conf1file, *conf2file, *matchndxfile, *outfile;
406 char title1[STRLEN],title2[STRLEN],*name1,*name2;
407 t_topology *top1,*top2;
409 t_atoms *atoms1,*atoms2;
412 real *w_rls,mass,totmass;
413 rvec *x1,*v1,*x2,*v2,*fit_x;
421 /* center of mass calculation */
425 /* variables for fit */
426 char *groupnames1,*groupnames2;
428 atom_id *index1,*index2;
429 real rms,msd,minmsd,maxmsd;
433 CopyRight(stderr,argv[0]);
434 parse_common_args(&argc,argv,PCA_BE_NICE | PCA_CAN_VIEW,
435 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
436 matchndxfile = opt2fn_null("-no",NFILE,fnm);
437 conf1file = ftp2fn(efTPS,NFILE,fnm);
438 conf2file = ftp2fn(efSTX,NFILE,fnm);
440 /* reading reference structure from first structure file */
441 fprintf(stderr,"\nReading first structure file\n");
443 read_tps_conf(conf1file,title1,top1,&ePBC1,&x1,&v1,box1,TRUE);
444 atoms1 = &(top1->atoms);
445 fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
446 title1,atoms1->nr,atoms1->nres);
449 rm_gropbc(atoms1,x1,box1);
451 fprintf(stderr,"Select group from first structure\n");
452 get_index(atoms1,opt2fn_null("-n1",NFILE,fnm),
453 1,&isize1,&index1,&groupnames1);
456 if (bFit && (isize1 < 3))
457 gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
459 /* reading second structure file */
460 fprintf(stderr,"\nReading second structure file\n");
462 read_tps_conf(conf2file,title2,top2,&ePBC2,&x2,&v2,box2,TRUE);
463 atoms2 = &(top2->atoms);
464 fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
465 title2,atoms2->nr,atoms2->nres);
468 rm_gropbc(atoms2,x2,box2);
470 fprintf(stderr,"Select group from second structure\n");
471 get_index(atoms2,opt2fn_null("-n2",NFILE,fnm),
472 1,&isize2,&index2,&groupnames2);
475 find_matching_names(&isize1, index1, atoms1, &isize2, index2, atoms2);
477 fp = ffopen(matchndxfile,"w");
478 fprintf(fp, "; Matching atoms between %s from %s and %s from %s\n",
479 groupnames1, conf1file, groupnames2, conf2file);
480 fprintf(fp, "[ Match_%s_%s ]\n", conf1file, groupnames1);
481 for(i=0; i<isize1; i++)
482 fprintf(fp, "%4u%s", index1[i]+1, (i%15==14 || i==isize1-1)?"\n":" ");
483 fprintf(fp, "[ Match_%s_%s ]\n", conf2file, groupnames2);
484 for(i=0; i<isize2; i++)
485 fprintf(fp, "%4u%s", index2[i]+1, (i%15==14 || i==isize2-1)?"\n":" ");
489 /* check isizes, must be equal */
490 if ( isize2 != isize1 )
491 gmx_fatal(FARGS,"You selected groups with differen number of atoms.\n");
493 for(i=0; i<isize1; i++) {
494 name1=*atoms1->atomname[index1[i]];
495 name2=*atoms2->atomname[index2[i]];
496 if (strcmp(name1,name2)) {
499 "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
500 i+1,index1[i]+1,name1,index2[i]+1,name2);
504 atoms1->atom[index1[i]].m = 1;
505 atoms2->atom[index2[i]].m = 1;
509 fprintf(stderr,"%d atomname%s did not match\n",warn,(warn==1) ? "":"s");
512 /* calculate and remove center of mass of structures */
513 calc_rm_cm(isize1, index1, atoms1, x1, xcm1);
514 calc_rm_cm(isize2, index2, atoms2, x2, xcm2);
516 snew(w_rls,atoms2->nr);
517 snew(fit_x,atoms2->nr);
518 for(at=0; (at<isize1); at++) {
519 w_rls[index2[at]] = atoms1->atom[index1[at]].m;
520 copy_rvec(x1[index1[at]],fit_x[index2[at]]);
523 /* do the least squares fit to the reference structure */
524 do_fit(atoms2->nr,w_rls,fit_x,x2);
536 /* calculate the rms deviation */
542 for(at=0; at<isize1; at++) {
543 mass = atoms1->atom[index1[at]].m;
544 for(m=0; m<DIM; m++) {
545 msd = sqr(x1[index1[at]][m] - x2[index2[at]][m]);
549 maxmsd = max(maxmsd, msds[at]);
550 minmsd = min(minmsd, msds[at]);
553 rms = sqrt(rms/totmass);
555 printf("Root mean square deviation after lsq fit = %g nm\n",rms);
557 printf("Atomic MSD's range from %g to %g nm^2\n",minmsd, maxmsd);
560 /* reset coordinates of reference and fitted structure */
561 for(i=0; i<atoms1->nr; i++)
564 for(i=0; i<atoms2->nr; i++)
569 outfile = ftp2fn(efSTO,NFILE,fnm);
570 switch (fn2ftp(outfile)) {
574 if (bBfac || bLabel) {
575 srenew(atoms1->pdbinfo, atoms1->nr);
576 srenew(atoms1->atom, atoms1->nr); /* Why renew atom? */
578 /* Avoid segfaults when writing the pdb-file */
579 for (i=0; i<atoms1->nr; i++) {
580 atoms1->pdbinfo[i].type = eptAtom;
581 atoms1->pdbinfo[i].occup = 1.00;
582 atoms1->pdbinfo[i].bAnisotropic = FALSE;
584 atoms1->pdbinfo[i].bfac = 0;
586 atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';
589 for(i=0; i<isize1; i++) {
590 /* atoms1->pdbinfo[index1[i]].type = eptAtom; */
591 /* atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */
593 atoms1->pdbinfo[index1[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
595 /* atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */
597 srenew(atoms2->pdbinfo, atoms2->nr);
598 srenew(atoms2->atom, atoms2->nr); /* Why renew atom? */
600 for (i=0; i<atoms2->nr; i++) {
601 atoms2->pdbinfo[i].type = eptAtom;
602 atoms2->pdbinfo[i].occup = 1.00;
603 atoms2->pdbinfo[i].bAnisotropic = FALSE;
605 atoms2->pdbinfo[i].bfac = 0;
607 atoms2->resinfo[atoms1->atom[i].resind].chainid = 'B';
610 for(i=0; i<isize2; i++) {
611 /* atoms2->pdbinfo[index2[i]].type = eptAtom; */
612 /* atoms2->pdbinfo[index2[i]].bAnisotropic = FALSE; */
614 atoms2->pdbinfo[index2[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
616 /* atoms2->resinfo[atoms2->atom[index2[i]].resind].chain = 'B'; */
619 fp=ffopen(outfile,"w");
621 write_pdbfile(fp,title1,atoms1,x1,ePBC1,box1,' ',1,NULL,TRUE);
622 write_pdbfile(fp,title2,atoms2,x2,ePBC2,box2,' ',bOne ? -1 : 2,NULL,TRUE);
627 fprintf(stderr,"WARNING: cannot write B-factor values to gro file\n");
628 fp=ffopen(outfile,"w");
630 write_hconf_p(fp,title1,atoms1,3,x1,v1,box1);
631 write_hconf_p(fp,title2,atoms2,3,x2,v2,box2);
636 fprintf(stderr,"WARNING: cannot write B-factor values to %s file\n",
637 ftp2ext(fn2ftp(outfile)));
640 "WARNING: cannot write the reference structure to %s file\n",
641 ftp2ext(fn2ftp(outfile)));
642 write_sto_conf(outfile,title2,atoms2,x2,v2,ePBC2,box2);
646 view_all(oenv,NFILE, fnm);