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 * GROningen Mixture of Alchemy and Childrens' Stories
42 #include "gmx_random.h"
53 void check_binwidth(real binwidth
) {
54 real smallest_bin
=0.1;
55 if (binwidth
<smallest_bin
)
56 gmx_fatal(FARGS
,"Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
59 void check_mcover(real mcover
) {
61 gmx_fatal(FARGS
,"mcover should be -1 or (0,1]");
62 } else if ((mcover
<0)&(mcover
!= -1)) {
63 gmx_fatal(FARGS
,"mcover should be -1 or (0,1]");
69 void normalize_probability(int n
,double *a
){
72 for (i
=0;i
<n
;i
++) norm
+=a
[i
];
73 for (i
=0;i
<n
;i
++) a
[i
]/=norm
;
76 gmx_nentron_atomic_structurefactors_t
*gmx_neutronstructurefactors_init(const char *datfn
) {
77 /* read nsfactor.dat */
85 gmx_nentron_atomic_structurefactors_t
*gnsf
;
89 /* allocate memory for structure */
91 snew(gnsf
->atomnm
,nralloc
);
92 snew(gnsf
->p
,nralloc
);
93 snew(gnsf
->n
,nralloc
);
94 snew(gnsf
->slength
,nralloc
);
96 gnsf
->nratoms
=line_no
;
98 while(get_a_line(fp
,line
,STRLEN
)) {
100 if (sscanf(line
,"%s %d %d %lf",atomnm
,&p
,&n
,&slength
) == 4) {
101 gnsf
->atomnm
[i
]=strdup(atomnm
);
104 gnsf
->slength
[i
]=slength
;
106 gnsf
->nratoms
=line_no
;
107 if (line_no
==nralloc
){
109 srenew(gnsf
->atomnm
,nralloc
);
110 srenew(gnsf
->p
,nralloc
);
111 srenew(gnsf
->n
,nralloc
);
112 srenew(gnsf
->slength
,nralloc
);
115 fprintf(stderr
,"WARNING: Error in file %s at line %d ignored\n",
118 srenew(gnsf
->atomnm
,gnsf
->nratoms
);
119 srenew(gnsf
->p
,gnsf
->nratoms
);
120 srenew(gnsf
->n
,gnsf
->nratoms
);
121 srenew(gnsf
->slength
,gnsf
->nratoms
);
125 return (gmx_nentron_atomic_structurefactors_t
*) gnsf
;
128 gmx_sans_t
*gmx_sans_init (t_topology
*top
, gmx_nentron_atomic_structurefactors_t
*gnsf
) {
129 gmx_sans_t
*gsans
=NULL
;
131 /* Try to assing scattering length from nsfactor.dat */
133 snew(gsans
->slength
,top
->atoms
.nr
);
134 /* copy topology data */
136 for(i
=0;i
<top
->atoms
.nr
;i
++) {
137 for(j
=0;j
<gnsf
->nratoms
;j
++) {
138 if(top
->atoms
.atom
[i
].atomnumber
== gnsf
->p
[j
]) {
139 /* we need special case for H and D */
140 if(top
->atoms
.atom
[i
].atomnumber
== 1) {
141 if(top
->atoms
.atom
[i
].m
== 1.008000) {
142 gsans
->slength
[i
] = gnsf
->slength
[0];
144 gsans
->slength
[i
] = gnsf
->slength
[1];
146 gsans
->slength
[i
] = gnsf
->slength
[j
];
151 return (gmx_sans_t
*) gsans
;
154 gmx_radial_distribution_histogram_t
*calc_radial_distribution_histogram (
164 gmx_radial_distribution_histogram_t
*pr
=NULL
;
172 gmx_rng_t
*trng
=NULL
;
174 gmx_large_int_t mc
=0,max
;
177 /* allocate memory for pr */
179 /* set some fields */
180 pr
->binwidth
=binwidth
;
183 * create max dist rvec
184 * dist = box[xx] + box[yy] + box[zz]
186 rvec_add(box
[XX
],box
[YY
],dist
);
187 rvec_add(box
[ZZ
],dist
,dist
);
191 pr
->grn
=(int)floor(rmax
/pr
->binwidth
)+1;
192 rmax
=pr
->grn
*pr
->binwidth
;
194 snew(pr
->gr
,pr
->grn
);
197 /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
199 max
=(gmx_large_int_t
)floor(0.5*0.01*isize
*(isize
-1));
201 max
=(gmx_large_int_t
)floor(0.5*mcover
*isize
*(isize
-1));
203 rng
=gmx_rng_init(seed
);
205 nthreads
= omp_get_max_threads();
208 for(i
=0;i
<nthreads
;i
++){
209 snew(tgr
[i
],pr
->grn
);
210 trng
[i
]=gmx_rng_init(gmx_rng_uniform_uint32(rng
));
212 #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
214 tid
= omp_get_thread_num();
215 /* now starting parallel threads */
217 for(mc
=0;mc
<max
;mc
++) {
218 i
=(int)floor(gmx_rng_uniform_real(trng
[tid
])*isize
);
219 j
=(int)floor(gmx_rng_uniform_real(trng
[tid
])*isize
);
221 tgr
[tid
][(int)floor(sqrt(distance2(x
[index
[i
]],x
[index
[j
]]))/binwidth
)]+=gsans
->slength
[index
[i
]]*gsans
->slength
[index
[j
]];
225 /* collecting data from threads */
226 for(i
=0;i
<pr
->grn
;i
++) {
227 for(j
=0;j
<nthreads
;j
++) {
228 pr
->gr
[i
] += tgr
[j
][i
];
231 /* freeing memory for tgr and destroying trng */
232 for(i
=0;i
<nthreads
;i
++) {
234 gmx_rng_destroy(trng
[i
]);
239 for(mc
=0;mc
<max
;mc
++) {
240 i
=(int)floor(gmx_rng_uniform_real(rng
)*isize
);
241 j
=(int)floor(gmx_rng_uniform_real(rng
)*isize
);
243 pr
->gr
[(int)floor(sqrt(distance2(x
[index
[i
]],x
[index
[j
]]))/binwidth
)]+=gsans
->slength
[index
[i
]]*gsans
->slength
[index
[j
]];
246 gmx_rng_destroy(rng
);
249 nthreads
= omp_get_max_threads();
250 /* Allocating memory for tgr arrays */
252 for(i
=0;i
<nthreads
;i
++) {
253 snew(tgr
[i
],pr
->grn
);
255 #pragma omp parallel shared(tgr) private(tid,i,j)
257 tid
= omp_get_thread_num();
258 /* starting parallel threads */
260 for(i
=0;i
<isize
;i
++) {
262 tgr
[tid
][(int)floor(sqrt(distance2(x
[index
[i
]],x
[index
[j
]]))/binwidth
)]+=gsans
->slength
[index
[i
]]*gsans
->slength
[index
[j
]];
266 /* collecating data for pr->gr */
267 for(i
=0;i
<pr
->grn
;i
++) {
268 for(j
=0;j
<nthreads
;j
++) {
269 pr
->gr
[i
] += tgr
[j
][i
];
272 /* freeing memory for tgr */
273 for(i
=0;i
<nthreads
;i
++) {
278 for(i
=0;i
<isize
;i
++) {
280 pr
->gr
[(int)floor(sqrt(distance2(x
[index
[i
]],x
[index
[j
]]))/binwidth
)]+=gsans
->slength
[index
[i
]]*gsans
->slength
[index
[j
]];
287 normalize_probability(pr
->grn
,pr
->gr
);
289 for(i
=0;i
<pr
->grn
;i
++)
290 pr
->r
[i
]=(pr
->binwidth
*i
+pr
->binwidth
*0.5);
292 return (gmx_radial_distribution_histogram_t
*) pr
;
295 gmx_static_structurefator_t
*convert_histogram_to_intensity_curve (gmx_radial_distribution_histogram_t
*pr
, double start_q
, double end_q
, double q_step
) {
296 gmx_static_structurefator_t
*sq
=NULL
;
300 sq
->qn
=(int)floor((end_q
-start_q
)/q_step
);
303 for(i
=0;i
<sq
->qn
;i
++)
304 sq
->q
[i
]=start_q
+i
*q_step
;
308 for(i
=1;i
<sq
->qn
;i
++) {
309 for(j
=0;j
<pr
->grn
;j
++)
310 sq
->s
[i
]+=(pr
->gr
[j
]/pr
->r
[j
])*sin(sq
->q
[i
]*pr
->r
[j
]);
311 sq
->s
[i
] /= sq
->q
[i
];
314 for(i
=0;i
<sq
->qn
;i
++) {
315 for(j
=0;j
<pr
->grn
;j
++)
316 sq
->s
[i
]+=(pr
->gr
[j
]/pr
->r
[j
])*sin(sq
->q
[i
]*pr
->r
[j
]);
317 sq
->s
[i
] /= sq
->q
[i
];
321 return (gmx_static_structurefator_t
*) sq
;