4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
56 int read_eigval(char *fn
,int nmax
,real eigval
[])
67 while (fgets(line
,STRLEN
-1,fp
) && !bEndOfSet
) {
68 bEndOfSet
= (line
[0] == '&');
69 if ((line
[0] != '#') && (line
[0] != '@') && !bEndOfSet
) {
70 if ((sscanf(line
,"%d %lf",&num
,&dbl
) != 2) || ((num
< 1) || (num
> nmax
)))
71 fprintf(stderr
,"Invalid line in %s: '%s'\n",fn
,line
);
83 int gmx_nmens(int argc
,char *argv
[])
85 static char *desc
[] = {
86 "[TT]g_nmens[tt] generates an ensemble around an average structure",
87 "in a subspace which is defined by a set of normal modes (eigenvectors).",
88 "The eigenvectors are assumed to be mass-weighted.",
89 "The position along each eigenvector is randomly taken from a Gaussian",
90 "distribution with variance kT/eigenvalue.[PAR]",
91 "By default the starting eigenvector is set to 7, since the first six",
92 "normal modes are the translational and rotational degrees of freedom."
94 static int nstruct
=100,first
=7,last
=-1,seed
=-1;
95 static real temp
=300.0;
97 { "-temp", FALSE
, etREAL
, {&temp
},
98 "Temperature in Kelvin" },
99 { "-seed", FALSE
, etINT
, {&seed
},
100 "Random seed, -1 generates a seed from time and pid" },
101 { "-num", FALSE
, etINT
, {&nstruct
},
102 "Number of structures to generate" },
103 { "-first", FALSE
, etINT
, {&first
},
104 "First eigenvector to use (-1 is select)" },
105 { "-last", FALSE
, etINT
, {&last
},
106 "Last eigenvector to use (-1 is till the last)" }
108 #define NPA asize(pa)
114 rvec
*xtop
,*xref
,*xav
,*xout1
,*xout2
;
116 int nvec
,*eignr
=NULL
;
119 real
*eigval
,totmass
,*invsqrtm
,t
,disp
;
121 char *grpname
,*indexfile
,title
[STRLEN
];
123 int nout
,*iout
,noutvec
,*outvec
;
125 real rfac
,invfr
,rhalf
,jr
;
127 const unsigned long im
= 0xffff;
128 const unsigned long ia
= 1093;
129 const unsigned long ic
= 18257;
132 { efTRN
, "-v", "eigenvec", ffREAD
},
133 { efXVG
, "-e", "eigenval", ffREAD
},
134 { efTPS
, NULL
, NULL
, ffREAD
},
135 { efNDX
, NULL
, NULL
, ffOPTRD
},
136 { efTRX
, "-o", "ensemble", ffWRITE
}
138 #define NFILE asize(fnm)
140 CopyRight(stderr
,argv
[0]);
141 parse_common_args(&argc
,argv
,PCA_BE_NICE
,
142 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,0,NULL
);
144 indexfile
=ftp2fn_null(efNDX
,NFILE
,fnm
);
146 read_eigenvectors(opt2fn("-v",NFILE
,fnm
),&natoms
,&bFit
,
147 &xref
,&bDMR
,&xav
,&bDMA
,&nvec
,&eignr
,&eigvec
);
149 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&xtop
,NULL
,box
,bDMA
);
152 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms
);
153 get_index(atoms
,indexfile
,1,&i
,&index
,&grpname
);
155 fatal_error(0,"you selected a group with %d elements instead of %d",
159 snew(eigval
,DIM
*natoms
);
160 neigval
=read_eigval(ftp2fn(efXVG
,NFILE
,fnm
),DIM
*natoms
,eigval
);
161 fprintf(stderr
,"Read %d eigenvalues\n",neigval
);
163 snew(invsqrtm
,natoms
);
165 for(i
=0; (i
<natoms
); i
++)
166 invsqrtm
[i
] = invsqrt(atoms
->atom
[index
[i
]].m
);
168 for(i
=0; (i
<natoms
); i
++)
175 /* make an index from first to last */
178 for(i
=0; i
<nout
; i
++)
181 printf("Select eigenvectors for output, end your selection with 0\n");
187 scanf("%d",&iout
[nout
]);
189 } while (iout
[nout
]>=0);
192 /* make an index of the eigenvectors which are present */
195 for(i
=0; i
<nout
; i
++) {
197 while ((j
<nvec
) && (eignr
[j
]!=iout
[i
]))
199 if ((j
<nvec
) && (eignr
[j
]==iout
[i
])) {
201 iout
[noutvec
] = iout
[i
];
205 fprintf(stderr
,"%d eigenvectors selected for output\n",noutvec
);
209 fprintf(stderr
,"Using seed %d and a temperature of %g K\n",seed
,temp
);
212 snew(xout2
,atoms
->nr
);
213 out
=open_trx(ftp2fn(efTRX
,NFILE
,fnm
),"w");
214 jran
= (unsigned long)((real
)im
*rando(&seed
));
215 for(s
=0; s
<nstruct
; s
++) {
216 for(i
=0; i
<natoms
; i
++)
217 copy_rvec(xav
[i
],xout1
[i
]);
218 for(j
=0; j
<noutvec
; j
++) {
220 /* (r-0.5) n times: var_n = n * var_1 = n/12
221 n=4: var_n = 1/3, so multiply with 3 */
223 rfac
= sqrt(3.0 * BOLTZ
*temp
/eigval
[iout
[j
]]);
225 rfac
= rfac
/(real
)im
;
227 jran
= (jran
*ia
+ic
) & im
;
229 jran
= (jran
*ia
+ic
) & im
;
231 jran
= (jran
*ia
+ic
) & im
;
233 jran
= (jran
*ia
+ic
) & im
;
235 disp
= rfac
* jr
- rhalf
;
237 for(i
=0; i
<natoms
; i
++)
239 xout1
[i
][d
] += disp
*eigvec
[v
][i
][d
]*invsqrtm
[i
];
241 for(i
=0; i
<natoms
; i
++)
242 copy_rvec(xout1
[i
],xout2
[index
[i
]]);
244 write_trx(out
,natoms
,index
,atoms
,0,t
,box
,xout2
,NULL
);
245 fprintf(stderr
,"\rGenerated %d structures",s
+1);
247 fprintf(stderr
,"\n");