3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
48 #include "gmx_fatal.h"
51 static real
dointerp(int n
,rvec x1
[],rvec x2
[],rvec xx
[],
52 int I
,int N
,real first
,real last
)
57 fac
= first
+ (I
*(last
-first
))/(N
-1);
61 for(j
=0; (j
<DIM
); j
++)
62 xx
[i
][j
] = fac0
*x1
[i
][j
] + fac1
*x2
[i
][j
];
67 int gmx_morph(int argc
,char *argv
[])
69 const char *desc
[] = {
70 "g_morph does a linear interpolation of conformations in order to",
71 "create intermediates. Of course these are completely unphysical, but",
72 "that you may try to justify yourself. Output is in the form of a ",
73 "generic trajectory. The number of intermediates can be controlled with",
74 "the -ninterm flag. The first and last flag correspond to the way of",
75 "interpolating: 0 corresponds to input structure 1 while",
76 "1 corresponds to input structure 2.",
77 "If you specify first < 0 or last > 1 extrapolation will be",
78 "on the path from input structure x1 to x2. In general the coordinates",
79 "of the intermediate x(i) out of N total intermidates correspond to:[PAR]",
80 "x(i) = x1 + (first+(i/(N-1))*(last-first))*(x2-x1)[PAR]",
81 "Finally the RMSD with respect to both input structures can be computed",
82 "if explicitly selected (-or option). In that case an index file may be",
83 "read to select what group RMS is computed from."
86 { efSTX
, "-f1", "conf1", ffREAD
},
87 { efSTX
, "-f2", "conf2", ffREAD
},
88 { efTRX
, "-o", "interm", ffWRITE
},
89 { efXVG
, "-or", "rms-interm", ffOPTWR
},
90 { efNDX
, "-n", "index", ffOPTRD
}
92 #define NFILE asize(fnm)
93 static int ninterm
= 11;
94 static real first
= 0.0;
95 static real last
= 1.0;
96 static gmx_bool bFit
= TRUE
;
98 { "-ninterm", FALSE
, etINT
, {&ninterm
},
99 "Number of intermediates" },
100 { "-first", FALSE
, etREAL
, {&first
},
101 "Corresponds to first generated structure (0 is input x0, see above)" },
102 { "-last", FALSE
, etREAL
, {&last
},
103 "Corresponds to last generated structure (1 is input x1, see above)" },
104 { "-fit", FALSE
, etBOOL
, {&bFit
},
105 "Do a least squares fit of the second to the first structure before interpolating" }
107 const char *leg
[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
109 int i
,isize
,is_lsq
,nat1
,nat2
;
111 atom_id
*index
,*index_lsq
,*index_all
,*dummy
;
115 real rms1
,rms2
,fac
,*mass
;
116 char title
[STRLEN
],*grpname
;
120 CopyRight(stderr
,argv
[0]);
121 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
,
122 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,
124 get_stx_coordnum (opt2fn("-f1",NFILE
,fnm
),&nat1
);
125 get_stx_coordnum (opt2fn("-f2",NFILE
,fnm
),&nat2
);
127 gmx_fatal(FARGS
,"Number of atoms in first structure is %d, in second %d",
130 init_t_atoms(&atoms
,nat1
,TRUE
);
136 read_stx_conf(opt2fn("-f1",NFILE
,fnm
),title
,&atoms
,x1
,v
,NULL
,box
);
137 read_stx_conf(opt2fn("-f2",NFILE
,fnm
),title
,&atoms
,x2
,v
,NULL
,box
);
140 snew(index_all
,nat1
);
141 for(i
=0; (i
<nat1
); i
++) {
146 printf("Select group for LSQ superposition:\n");
147 get_index(&atoms
,opt2fn_null("-n",NFILE
,fnm
),1,&is_lsq
,&index_lsq
,
149 reset_x(is_lsq
,index_lsq
,nat1
,index_all
,x1
,mass
);
150 reset_x(is_lsq
,index_lsq
,nat1
,index_all
,x2
,mass
);
151 do_fit(nat1
,mass
,x1
,x2
);
154 bRMS
= opt2bSet("-or",NFILE
,fnm
);
156 fp
= xvgropen(opt2fn("-or",NFILE
,fnm
),"RMSD","Conf","(nm)",oenv
);
157 xvgr_legend(fp
,asize(leg
),leg
,oenv
);
158 printf("Select group for RMSD calculation:\n");
159 get_index(&atoms
,opt2fn_null("-n",NFILE
,fnm
),1,&isize
,&index
,&grpname
);
160 printf("You selected group %s, containing %d atoms\n",grpname
,isize
);
161 rms1
= rmsdev_ind(isize
,index
,mass
,x1
,x2
);
162 fprintf(stderr
,"RMSD between input conformations is %g nm\n",rms1
);
166 for(i
=0; (i
<nat1
); i
++)
168 status
= open_trx(ftp2fn(efTRX
,NFILE
,fnm
),"w");
170 for(i
=0; (i
<ninterm
); i
++) {
171 fac
= dointerp(nat1
,x1
,x2
,xx
,i
,ninterm
,first
,last
);
172 write_trx(status
,nat1
,dummy
,&atoms
,i
,fac
,box
,xx
,NULL
,NULL
);
174 rms1
= rmsdev_ind(isize
,index
,mass
,x1
,xx
);
175 rms2
= rmsdev_ind(isize
,index
,mass
,x2
,xx
);
176 fprintf(fp
,"%10g %10g %10g\n",fac
,rms1
,rms2
);
184 do_view(oenv
,opt2fn("-or",NFILE
,fnm
),"-nxy");