4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
27 * GRowing Old MAkes el Chrono Sweat
29 static char *SRCID_g_covar_c
= "$Id$";
56 int main(int argc
,char *argv
[])
58 static char *desc
[] = {
59 "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
61 "All structures are fitted to the structure in the structure file.",
62 "When this is not a run input file periodicity will not be taken into",
63 "account. When the fit and analysis groups are identical and the analysis",
64 "is non mass-weighted, the fit will also be non mass-weighted.[PAR]",
65 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
66 "When the same atoms are used for the fit and the covariance analysis,",
67 "the reference structure is written first with t=-1.",
68 "The average structure is written with t=0, the eigenvectors",
69 "are written as frames with the eigenvector number as timestamp.",
70 "The eigenvectors can be analyzed with [TT]g_anaeig[tt]."
72 static bool bFit
=TRUE
,bM
=FALSE
;
75 { "-fit", FALSE
, etBOOL
, {&bFit
},
76 "Fit to a reference structure"},
77 { "-mwa", FALSE
, etBOOL
, {&bM
},
78 "Mass-weighted covariance analysis"},
79 { "-last", FALSE
, etINT
, {&end
},
80 "Last eigenvector to write away (-1 is till the last)" }
86 rvec
*x
,*xread
,*xref
,*xav
;
88 real t
,tstart
,tend
,*mat
,dev
,trace
,sum
,*eigval
,inv_nframes
;
89 real xj
,*sqrtm
,*w_rls
=NULL
;
91 int natoms
,nat
,ndim
,count
,nframes
;
93 char *fitfile
,*trxfile
,*ndxfile
;
94 char *eigvalfile
,*eigvecfile
,*averfile
,*logfile
;
95 char str
[STRLEN
],*fitname
,*ananame
;
96 int i
,j
,k
,l
,d
,dj
,nfit
;
97 atom_id
*index
,*all_at
,*ifit
;
98 bool bDiffMass1
,bDiffMass2
;
101 { efTRX
, "-f", NULL
, ffREAD
},
102 { efTPS
, NULL
, NULL
, ffREAD
},
103 { efNDX
, NULL
, NULL
, ffOPTRD
},
104 { efXVG
, NULL
, "eigenval", ffWRITE
},
105 { efTRN
, "-v", "eigenvec", ffWRITE
},
106 { efSTO
, "-av", "average.pdb", ffWRITE
},
107 { efLOG
, NULL
, "covar", ffWRITE
}
109 #define NFILE asize(fnm)
111 CopyRight(stderr
,argv
[0]);
112 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_SET_NPRI
,TRUE
,
113 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
);
117 fitfile
= ftp2fn(efTPS
,NFILE
,fnm
);
118 trxfile
= ftp2fn(efTRX
,NFILE
,fnm
);
119 ndxfile
= ftp2fn_null(efNDX
,NFILE
,fnm
);
120 eigvalfile
= ftp2fn(efXVG
,NFILE
,fnm
);
121 eigvecfile
= ftp2fn(efTRN
,NFILE
,fnm
);
122 averfile
= ftp2fn(efSTO
,NFILE
,fnm
);
123 logfile
= ftp2fn(efLOG
,NFILE
,fnm
);
125 read_tps_conf(fitfile
,str
,&top
,&xref
,NULL
,box
,TRUE
);
129 printf("\nChoose a group for the least squares fit\n");
130 get_index(atoms
,ndxfile
,1,&nfit
,&ifit
,&fitname
);
132 fatal_error(0,"Need >= 3 points to fit!\n");
135 printf("\nChoose a group for the covariance analysis\n");
136 get_index(atoms
,ndxfile
,1,&natoms
,&index
,&ananame
);
140 snew(w_rls
,atoms
->nr
);
141 for(i
=0; (i
<nfit
); i
++) {
142 w_rls
[ifit
[i
]]=atoms
->atom
[ifit
[i
]].m
;
144 bDiffMass1
= bDiffMass1
|| (w_rls
[ifit
[i
]]!=w_rls
[ifit
[i
-1]]);
149 for(i
=0; (i
<natoms
); i
++)
151 sqrtm
[i
]=sqrt(atoms
->atom
[index
[i
]].m
);
153 bDiffMass2
= bDiffMass2
|| (sqrtm
[i
]!=sqrtm
[i
-1]);
158 if (bFit
&& bDiffMass1
&& !bDiffMass2
) {
159 bDiffMass1
= natoms
!= nfit
;
161 for (i
=0; (i
<natoms
) && !bDiffMass1
; i
++)
162 bDiffMass1
= index
[i
] != ifit
[i
];
165 "Note: the fit and analysis group are identical,\n"
166 " while the fit is mass weighted and the analysis is not.\n"
167 " Making the fit non mass weighted.\n\n");
168 for(i
=0; (i
<nfit
); i
++)
173 snew(all_at
,atoms
->nr
);
174 for(i
=0; (i
<atoms
->nr
); i
++)
177 /* Prepare reference frame */
178 rm_pbc(&(top
.idef
),atoms
->nr
,box
,xref
,xref
);
180 reset_x(nfit
,ifit
,atoms
->nr
,all_at
,xref
,w_rls
);
187 fprintf(stderr
,"Constructing covariance matrix (%dx%d)...\n",ndim
,ndim
);
190 nat
=read_first_x(&status
,trxfile
,&t
,&xread
,box
);
195 /* calculate x: a fitted struture of the selected atoms */
196 rm_pbc(&(top
.idef
),nat
,box
,xread
,xread
);
198 reset_x(nfit
,ifit
,nat
,all_at
,xread
,w_rls
);
199 do_fit(nat
,w_rls
,xref
,xread
);
201 for (i
=0; i
<natoms
; i
++)
202 copy_rvec(xread
[index
[i
]],x
[i
]);
204 for (j
=0; j
<natoms
; j
++) {
205 /* calculate average structure */
206 rvec_inc(xav
[j
],x
[j
]);
207 /* calculate cross product matrix */
208 for (dj
=0; dj
<DIM
; dj
++) {
211 for (i
=j
; i
<natoms
; i
++) {
214 mat
[l
+d
] += x
[i
][d
]*xj
;
218 } while (read_next_x(status
,&t
,nat
,xread
,box
));
221 /* calculate the mass-weighted covariance matrix */
222 inv_nframes
=1.0/nframes
;
223 for (i
=0; i
<natoms
; i
++) {
224 svmul(inv_nframes
,xav
[i
],xav
[i
]);
225 copy_rvec(xav
[i
],xread
[index
[i
]]);
228 write_sto_conf_indexed(opt2fn("-av",NFILE
,fnm
),"Average structure",
229 atoms
,xread
,NULL
,NULL
,natoms
,index
);
231 for (j
=0; j
<natoms
; j
++)
232 for (dj
=0; dj
<DIM
; dj
++)
233 for (i
=j
; i
<natoms
; i
++) {
234 k
=ndim
*(DIM
*j
+dj
)+DIM
*i
;
235 for (d
=0; d
<DIM
; d
++) {
236 mat
[k
+d
]=(mat
[k
+d
]*inv_nframes
-xav
[i
][d
]*xav
[j
][dj
])
240 for (j
=0; j
<ndim
; j
++)
241 for (i
=j
; i
<ndim
; i
++)
242 mat
[ndim
*i
+j
]=mat
[ndim
*j
+i
];
245 for(i
=0; i
<ndim
; i
++)
246 trace
+=mat
[i
*ndim
+i
];
247 fprintf(stderr
,"\nTrace of the covariance matrix: %g (%snm^2)\n",
248 trace
,bM
? "amu " : "");
251 fprintf(stderr
,"Dumping the covariance matrix to covmat.dat\n");
252 out
= ffopen("covmat.dat","w");
253 for (j
=0; j
<ndim
; j
++) {
254 for (i
=0; i
<ndim
; i
++)
255 fprintf(out
," %g",mat
[ndim
*j
+i
]);
260 /* call diagonalization routine */
262 fprintf(stderr
,"\nDiagonalizing...\n");
266 ql77 (ndim
,mat
,eigval
);
268 /* now write the output */
271 for(i
=0; i
<ndim
; i
++)
273 fprintf(stderr
,"\nSum of the eigenvalues: %g (%snm^2)\n",
274 sum
,bM
? "amu " : "");
275 if (fabs(trace
-sum
)>0.01*trace
)
276 fprintf(stderr
,"\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
278 fprintf(stderr
,"\nWriting eigenvalues to %s\n",eigvalfile
);
280 sprintf(str
,"(%snm\\S2\\N)",bM
? "amu " : "");
281 out
=xvgropen(eigvalfile
,
282 "Eigenvalues of the covariance matrix",
283 "Eigenvector index",str
);
284 for (i
=0; (i
<ndim
); i
++)
285 fprintf (out
,"%10d %g\n",i
+1,eigval
[ndim
-1-i
]);
289 if (nframes
-1 < ndim
)
295 /* misuse lambda: 0/1 mass weighted analysis no/yes */
297 WriteXref
= eWXR_YES
;
298 for(i
=0; i
<nfit
; i
++)
299 copy_rvec(xref
[ifit
[i
]],x
[i
]);
303 /* misuse lambda: -1 for no fit */
304 WriteXref
= eWXR_NOFIT
;
307 write_eigenvectors(eigvecfile
,natoms
,mat
,TRUE
,1,end
,
308 WriteXref
,x
,bDiffMass1
,xav
,bM
);
310 out
= ffopen(logfile
,"w");
313 fprintf(out
,"Covariance analysis log, written %s\n",
315 fprintf(out
,"Program: %s\n",argv
[0]);
317 fprintf(out
,"Working directory: %s\n\n",str
);
319 fprintf(out
,"Read %d frames from %s (time %g to %g)\n",nframes
,trxfile
,
322 fprintf(out
,"Read reference structure for fit from %s\n",fitfile
);
324 fprintf(out
,"Read index groups from %s\n",ndxfile
);
327 fprintf(out
,"Analysis group is '%s' (%d atoms)\n",ananame
,natoms
);
329 fprintf(out
,"Fit group is '%s' (%d atoms)\n",fitname
,nfit
);
331 fprintf(out
,"No fit was used\n");
332 fprintf(out
,"Analysis is %smass weighted\n", bDiffMass2
? "":"non-");
334 fprintf(out
,"Fit is %smass weighted\n", bDiffMass1
? "":"non-");
335 fprintf(out
,"Diagonalized the %dx%d covariance matrix\n",ndim
,ndim
);
336 fprintf(out
,"Trace of the covariance matrix before diagonalizing: %g\n",
338 fprintf(out
,"Trace of the covariance matrix after diagonalizing: %g\n\n",
341 fprintf(out
,"Wrote %d eigenvalues to %s\n",ndim
,eigvalfile
);
342 if (WriteXref
== eWXR_YES
)
343 fprintf(out
,"Wrote reference structure to %s\n",eigvecfile
);
344 fprintf(out
,"Wrote average structure to %s and %s\n",averfile
,eigvecfile
);
345 fprintf(out
,"Wrote eigenvectors %d to %d to %s\n",1,end
,eigvecfile
);
349 fprintf(stderr
,"Wrote the log to %s\n",logfile
);