?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed the execution order in OpenMM so that the generated trajectory frames correspon...
[gromacs.git] / src / contrib / mkice.c
blobf98249080ea8f974e659a0a872c74184edee0793
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.3.99_development_20071104
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-2006, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Groningen Machine for Chemical Simulation
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <stdio.h>
40 #include <math.h>
41 #include "typedefs.h"
42 #include "statutil.h"
43 #include "copyrite.h"
44 #include "gmx_fatal.h"
45 #include "pdbio.h"
46 #include "macros.h"
47 #include "smalloc.h"
48 #include "vec.h"
49 #include "pbc.h"
50 #include "physics.h"
51 #include "names.h"
52 #include "txtdump.h"
53 #include "trnio.h"
54 #include "symtab.h"
55 #include "strdb.h"
56 #include "atomprop.h"
57 #include "confio.h"
58
59 #define TET 109.47
60 #define DCONS 0.117265878
61
62 typedef struct {
63 int n,aa[4];
64 } t_bbb;
65
66 static char *watname[] = { "OW ", "HW1", "HW2", "DW", "SW" };
67 static char *diamname[] = { "C1", "H2", "H3", "H4", "H5", "H2", "H3", "H4", "H5" };
68 static real qspc[] = { -0.82, 0.41, 0.41 };
69 static real qyaw[] = { 1.24588, 0.62134, 0.62134, 0.0, -2.48856 };
70 static real spc_lj[3][6] = {
71 { 2.6171e-3, 2.6331e-6, 0, 0, 0, 0 },
72 { 0, 0, 0, 0, 0, 0 },
73 { 0, 0, 0, 0, 0, 0 }
74 };
75 #define CHH 2e-9
76 #define CHS 2e-9
77 #define COS 2e-6
78 static real yaw_lj[5][10] = {
79 { 0, 0, 0, 0, 0, 0, 0, 0, 0, COS },
80 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
81 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
82 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
83 { 0, COS, 0, CHS, 0, CHS, 0, 0, 2.6e-3, 0 }
84 };
85
86 void unitcell(rvec x[],rvec box,bool bYaw,real odist,real hdist)
87 {
88 #define cx 0.81649658
89 #define cy 0.47140452
90 #define cy2 0.94280904
91 #define cz 0.33333333
92
93 rvec xx[24] = {
94 { 0, 0, 0 }, /* O1 */
95 { 0, 0, 1 }, /* H relative to Oxygen */
96 { cx, cy, -cz },
97 { cx, cy, -cz }, /* O2 */
98 { 0, 0, -1 }, /* H relative to Oxygen */
99 { cx,-cy, +cz },
100 { cx, cy+cy2, 0 }, /* O3 */
101 { -cx, cy, -cz }, /* H relative to Oxygen */
102 { 0, -cy2, -cz },
103 { 0, 2*cy+cy2, -cz }, /* O4 */
104 {-cx,-cy, +cz }, /* H relative to Oxygen */
105 { 0 , cy2, +cz },
106 { 0, 0, 1 }, /* O5 */
107 {-cx, cy, +cz }, /* H relative to Oxygen */
108 { 0 , -cy2, +cz },
109 { cx, cy, 1+cz }, /* O6 */
110 { -cx, -cy, -cz }, /* H relative to Oxygen */
111 { 0, cy2, -cz },
112 { cx, cy+cy2, 1 }, /* O7 */
113 { 0, 0, -1 }, /* H relative to Oxygen */
114 { cx, cy, +cz },
115 { 0, 2*cy+cy2,1+cz }, /* O8 */
116 { 0, 0, 1 }, /* H relative to Oxygen */
117 { cx, -cy, -cz }
118 };
119 int i,iin,iout,j,m;
120 rvec tmp,t2,dip;
121
122 clear_rvec(dip);
123 for(i=0; (i<8); i++) {
124 iin = 3*i;
125 if (bYaw)
126 iout = 5*i;
127 else
128 iout = iin;
129 svmul(odist,xx[iin],x[iout]);
130 svmul(-0.82,x[iout],t2);
131 rvec_inc(dip,t2);
132 for(j=1; (j<=2); j++) {
133 svmul(hdist,xx[iin+j],tmp);
134 rvec_add(x[iout],tmp,x[iout+j]);
135 svmul(0.41,x[iout+j],t2);
136 rvec_inc(dip,t2);
137 }
138 if (bYaw)
139 for(m=0; (m<DIM); m++)
140 x[iout+3][m] = x[iout+4][m] =
141 (1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
142 }
143
144 box[XX] = 2*cx;
145 box[YY] = 2*(cy2+cy);
146 box[ZZ] = 2*(1+cz);
147 for(i=0; (i<DIM); i++)
148 box[i] *= odist;
149
150 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
151 printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
152 }
153
154 void random_h_coords(int natmol,int nmol,rvec x[],rvec box,
155 bool bYaw,real odist,real hdist)
156 {
157 #define cx 0.81649658
158 #define cy 0.47140452
159 #define cy2 0.94280904
160 #define cz 0.33333333
161
162 rvec xx[24] = {
163 { 0, 0, 0 }, /* O1 */
164 { 0, 0, 1 }, /* H relative to Oxygen */
165 { cx, cy, -cz },
166 { cx, cy, -cz }, /* O2 */
167 { 0, 0, -1 }, /* H relative to Oxygen */
168 { cx,-cy, +cz },
169 { cx, cy+cy2, 0 }, /* O3 */
170 { -cx, cy, -cz }, /* H relative to Oxygen */
171 { 0, -cy2, -cz },
172 { 0, 2*cy+cy2, -cz }, /* O4 */
173 {-cx,-cy, +cz }, /* H relative to Oxygen */
174 { 0 , cy2, +cz },
175 { 0, 0, 1 }, /* O5 */
176 {-cx, cy, +cz }, /* H relative to Oxygen */
177 { 0 , -cy2, +cz },
178 { cx, cy, 1+cz }, /* O6 */
179 { -cx, -cy, -cz }, /* H relative to Oxygen */
180 { 0, cy2, -cz },
181 { cx, cy+cy2, 1 }, /* O7 */
182 { 0, 0, -1 }, /* H relative to Oxygen */
183 { cx, cy, +cz },
184 { 0, 2*cy+cy2,1+cz }, /* O8 */
185 { 0, 0, 1 }, /* H relative to Oxygen */
186 { cx, -cy, -cz }
187 };
188 int i,iin,iout,j,m;
189 rvec tmp,t2,dip;
190
191 clear_rvec(dip);
192 for(i=0; (i<nmol); i++) {
193 iin = natmol*i;
194 iout = iin;
195 svmul(odist,x[iin],x[iout]);
196 svmul(-0.82,x[iout],t2);
197 rvec_inc(dip,t2);
198 for(j=1; (j<=2); j++) {
199 svmul(hdist,xx[3*(i % 8)+j],tmp);
200 rvec_add(x[iout],tmp,x[iout+j]);
201 svmul(0.41,x[iout+j],t2);
202 rvec_inc(dip,t2);
203 }
204 }
205
206 box[XX] = 2*cx;
207 box[YY] = 2*(cy2+cy);
208 box[ZZ] = 2*(1+cz);
209 for(i=0; (i<DIM); i++)
210 box[i] *= odist;
211
212 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
213 printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
214 }
215
216 void unitcell_d(rvec x[],rvec box,real odist)
217 {
218 rvec cc[8] = {
219 { 0, 0, 0 }, /* C1 */
220 { cx, cy, -cz }, /* C2 */
221 { cx, cy+cy2, 0 }, /* C3 */
222 { 0, 2*cy+cy2, -cz }, /* C4 */
223 { 0, 0, 1 }, /* C5 */
224 { cx, cy, 1+cz }, /* C6 */
225 { cx, cy+cy2, 1 }, /* C7 */
226 { 0, 2*cy+cy2,1+cz }, /* C8 */
227 };
228 rvec hh[4] = {
229 { 0, 0, 1 }, /* H relative to C */
230 { cx, cy, -cz },
231 { cx, -cy, -cz },
232 {-cy2, 0, -cz }
233 };
234 int i,iin,iout,j,m;
235 rvec tmp,t2,dip;
236
237 clear_rvec(dip);
238 for(i=0; (i<8); i++) {
239 iin = i;
240 iout = i;
241 svmul(odist,cc[iin],x[iout]);
242 }
243 box[XX] = 2*cx;
244 box[YY] = 2*(cy2+cy);
245 box[ZZ] = 2*(1+cz);
246 for(i=0; (i<DIM); i++)
247 box[i] *= odist;
248
249 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
250 }
251
252 static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
253 bool bPBC,matrix box)
254 {
255 real od2 = odist*odist+1e-5;
256 t_pbc pbc;
257 t_bbb *bbb;
258 int i,j;
259 rvec dx;
260
261 if (bPBC)
262 set_pbc(&pbc,box);
263 snew(bbb,natoms);
264 for(i=0; (i<natoms); i++) {
265 for(j=i+1; (j<natoms); j++) {
266 if (bPBC)
267 pbc_dx(&pbc,x[i],x[j],dx);
268 else
269 rvec_sub(x[i],x[j],dx);
270 if (iprod(dx,dx) <= od2) {
271 bbb[i].aa[bbb[i].n++] = j;
272 bbb[j].aa[bbb[j].n++] = i;
273 }
274 }
275 }
276 if (debug)
277 #define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
278 for(i=0; (i<natoms); i++)
279 fprintf(debug,"bonds from %3d: %d %d %d %d\n",
280 i,PRB(1),PRB(2),PRB(3),PRB(4));
281 #undef PRB
282 return bbb;
283 }
284
285 static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
286 bool bPBC,matrix box)
287 {
288
289 int i,ib,j,k,l,m,nrm=0;
290 t_bbb *bbb;
291 bool *bRemove;
292 rvec dx;
293
294 do {
295 nrm = 0;
296 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
297
298 for(i=0; (i<a->nr); i++) {
299 if (bbb[i].n < 2) {
300 for(k=0; (k<bbb[i].n); k++) {
301 ib = bbb[i].aa[k];
302 for(j=0; (j<bbb[ib].n); j++)
303 if (bbb[ib].aa[j] == i)
304 break;
305 if (j == bbb[ib].n)
306 gmx_fatal(FARGS,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
307 for( ; (j<bbb[ib].n-1); j++)
308 bbb[ib].aa[j] = bbb[ib].aa[j+1];
309 bbb[ib].n--;
310 nrm++;
311 }
312 bbb[i].n = 0;
313 }
314 }
315
316 for(i=j=0; (i<a->nr); i++) {
317 if (bbb[i].n >= 2) {
318 copy_rvec(x[i],x[j]);
319 j++;
320 }
321 }
322 fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
323 a->nr-j,a->nr);
324 a->nr = j;
325 sfree(bbb);
326 } while (nrm > 0);
327 /* Rename atoms */
328 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
329 for(i=0; (i<a->nr); i++) {
330 switch (bbb[i].n) {
331 case 4:
332 a->atomname[i] = put_symtab(symtab,"C");
333 break;
334 case 3:
335 a->atomname[i] = put_symtab(symtab,"CH1");
336 break;
337 case 2:
338 a->atomname[i] = put_symtab(symtab,"CH2");
339 break;
340 default:
341 gmx_fatal(FARGS,"This atom (%d) has %d bonds only",i,bbb[i].n);
342 }
343 }
344 sfree(bbb);
345 }
346
347 static real calc_ener(real c6,real c12,rvec dx,tensor vir)
348 {
349 real r,e,f;
350 int m,n;
351
352 r = norm(dx);
353 e = c12*pow(r,-12.0) - c6*pow(r,-6.0);
354 f = 12*c12*pow(r,-14.0) - 6*c6*pow(r,-8.0);
355 for(m=0; (m<DIM); m++)
356 for(n=0; (n<DIM); n++)
357 vir[m][n] -= 0.5*f*dx[m]*dx[n];
358
359 return e;
360 }
361
362 static int read_rel_coords(char *fn,rvec **xx,int natmol)
363 {
364 int i,nline;
365 double x,y,z;
366 char **strings=NULL;
367
368 nline = get_file(fn,&strings);
369 printf("Read %d lines from %s\n",nline,fn);
370 snew(*xx,nline*natmol);
371 for(i=0; (i<nline); i++) {
372 if (sscanf(strings[i],"%lf%lf%lf",&x,&y,&z) != 3)
373 gmx_fatal(FARGS,"Not enough arguments on line %d of file %s (should be 3)",i,fn);
374 (*xx)[natmol*i][XX] = x;
375 (*xx)[natmol*i][YY] = y;
376 (*xx)[natmol*i][ZZ] = z;
377 }
378 return natmol*nline;
379 }
380
381 void virial(FILE *fp,bool bFull,int nmol,rvec x[],matrix box,real rcut,
382 bool bYaw,real q[],bool bLJ)
383 {
384 int i,j,im,jm,natmol,ik,jk,m,ninter;
385 rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
386 real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
387 t_pbc pbc;
388
389 set_pbc(&pbc,box);
390 fprintf(fp,"%3s - %3s: %6s %6s %6s %6s %8s %8s %8s\n",
391 "ai","aj","dx","dy","dz","|d|","virx","viry","virz");
392 clear_rvec(ftot);
393 clear_rvec(vir);
394 ninter = 0;
395 vctot = 0;
396 v12tot = 0;
397 v6tot = 0;
398 natmol = bYaw ? 5 : 3;
399 snew(force,nmol*natmol);
400
401 for(i=0; (i<nmol); i++) {
402 im = natmol*i;
403 /* Center of geometry */
404 clear_rvec(xcmi);
405 for(ik=0; (ik<natmol); ik++)
406 rvec_inc(xcmi,x[im+ik]);
407 for(m=0; (m<DIM); m++)
408 xcmi[m] /= natmol;
409
410 for(j=i+1; (j<nmol); j++) {
411 jm = natmol*j;
412 /* Center of geometry */
413 clear_rvec(xcmj);
414 for(jk=0; (jk<natmol); jk++)
415 rvec_inc(xcmj,x[jm+jk]);
416 for(m=0; (m<DIM); m++)
417 xcmj[m] /= natmol;
418
419 /* First check COM-COM distance */
420 pbc_dx(&pbc,xcmi,xcmj,dx);
421 if (norm(dx) < rcut) {
422 ninter++;
423 /* Neirest neighbour molecules! */
424 clear_rvec(dvir);
425 for(ik=0; (ik<natmol); ik++) {
426 for(jk=0; (jk<natmol); jk++) {
427 pbc_dx(&pbc,x[im+ik],x[jm+jk],dx);
428 dx2 = iprod(dx,dx);
429 dx1 = sqrt(dx2);
430 vcoul = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
431 vctot += vcoul;
432
433 if (bLJ) {
434 if (bYaw) {
435 c6 = yaw_lj[ik][2*jk];
436 c12 = yaw_lj[ik][2*jk+1];
437 }
438 else {
439 c6 = spc_lj[ik][2*jk];
440 c12 = spc_lj[ik][2*jk+1];
441 }
442 dx6 = dx2*dx2*dx2;
443 v6 = c6/dx6;
444 v12 = c12/(dx6*dx6);
445 v6tot -= v6;
446 v12tot+= v12;
447 fscal = (vcoul+12*v12-6*v6)/dx2;
448 }
449 else
450 fscal = vcoul/dx2;
451 for(m=0; (m<DIM); m++) {
452 f[m] = dx[m]*fscal;
453 dvir[m] -= 0.5*dx[m]*f[m];
454 }
455 rvec_inc(force[ik+im],f);
456 rvec_dec(force[jk+jm],f);
457 /*if (bFull)
458 fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
459 " %8.3f %8.3f %8.3f\n",
460 watname[ik],im+ik,watname[jk],jm+jk,
461 dx[XX],dx[YY],dx[ZZ],norm(dx),
462 dvir[XX],dvir[YY],dvir[ZZ]);*/
463 }
464 }
465 if (bFull)
466 fprintf(fp,"%3s%4d-%3s%4d: "
467 " %8.3f %8.3f %8.3f\n",
468 "SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
469 rvec_inc(vir,dvir);
470 }
471 }
472 }
473 fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
474 ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
475 fprintf(fp,"Vcoul: %10.4e V12: %10.4e V6: %10.4e Vtot: %10.4e (kJ/mol)\n",
476 vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
477 pr_rvec(fp,0,"vir ",vir,DIM,TRUE);
478
479 for(m=0; (m<DIM); m++)
480 pres[m] = -2*PRESFAC/(det(box))*vir[m];
481 pr_rvec(fp,0,"pres",pres,DIM,TRUE);
482 pr_rvecs(fp,0,"force",force,natmol*nmol);
483 sfree(force);
484 }
485
486
487
488 int main(int argc,char *argv[])
489 {
490 static char *desc[] = {
491 "mkice generates an ice crystal in the Ih crystal form which is the",
492 "most stable form. The rectangular unitcell contains eight molecules",
493 "and all oxygens are tetrahedrally coordinated.[PAR]",
494 "If an input file is given it is interpreted as a series of oxygen",
495 "coordinates the distance between which can be scaled by the odist flag.",
496 "The program then adds hydrogens to the oxygens in random orientation",
497 "but with proper OH distances and HOH angle. This feature allows to",
498 "build water clusters based on oxygen coordinates only."
499 };
500 static int nx=1,ny=1,nz=1;
501 static bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE;
502 static bool bOrdered=TRUE,bDiamond=FALSE,bPBC=TRUE;
503 static real rcut=0.3,odist=0.274,hdist=0.09572;
504 t_pargs pa[] = {
505 { "-nx", FALSE, etINT, {&nx}, "nx" },
506 { "-ny", FALSE, etINT, {&ny}, "ny" },
507 { "-nz", FALSE, etINT, {&nz}, "nz" },
508 { "-yaw", FALSE, etBOOL, {&bYaw},
509 "Generate virtual sites and shell positions" },
510 { "-lj", FALSE, etBOOL, {&bLJ},
511 "Use LJ as well as coulomb for virial calculation" },
512 { "-rcut", FALSE,etREAL, {&rcut},"Cut-off for virial calculations" },
513 { "-full", FALSE,etBOOL, {&bFull},"Full virial output" },
514 { "-odist", FALSE, etREAL, {&odist}, "Distance (nm) between oxygens" },
515 { "-hdist", FALSE, etREAL, {&hdist}, "Bondlength (nm) for OH bond" },
516 { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
517 { "-pbc", FALSE, etBOOL, {&bPBC}, "Make a periodic diamond" },
518 { "-order", FALSE,etBOOL, {&bOrdered}, "Make a proton-ordered ice lattice" },
519 { "-series",FALSE, etBOOL, {&bSeries},
520 "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
521 };
522 t_filenm fnm[] = {
523 { efSTO, "-p", "ice", ffWRITE },
524 { efSTO, "-c", NULL, ffOPTRD },
525 { efDAT, "-f", NULL, ffOPTRD },
526 { efTRN, "-o", "ice", ffOPTWR }
527 };
528 #define NFILE asize(fnm)
529
530 FILE *fp;
531 char *fn,quote[256];
532 int i,j,k,n,nmax,m,natom,natmol;
533 t_atoms *pdba;
534 t_atoms atoms;
535 t_symtab symtab;
536 rvec box,tmp,*xx;
537 matrix boxje;
538
539 CopyRight(stdout,argv[0]);
540 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),
541 desc,0,NULL);
542 if (debug) {
543 fprintf(debug,"nx = %3d, ny = %3d, nz = %3d\n",nx,ny,nz);
544 fprintf(debug,"YAW = %3s, LJ = %3s, rcut = %g\n",yesno_names[bYaw],
545 yesno_names[bLJ],rcut);
546 }
547
548 if (bYaw)
549 natmol = 5;
550 else if (bDiamond)
551 natmol = 1;
552 else
553 natmol = 3;
554
555 if (opt2bSet("-f",NFILE,fnm)) {
556 natom = read_rel_coords(opt2fn("-f",NFILE,fnm),&xx,natmol);
557 nmax = natom;
558 }
559 else {
560 natom = natmol*8;
561 nmax = natom*nx*ny*nz;
562 snew(xx,nmax);
563 }
564 snew(pdba,1);
565 init_t_atoms(pdba,nmax,TRUE);
566 pdba->nr = nmax;
567 open_symtab(&symtab);
568 for(n=0; (n<nmax); n++) {
569 pdba->pdbinfo[n].type = epdbATOM;
570 pdba->pdbinfo[n].atomnr = 1+n;
571 pdba->atom[n].resnr = 1+(n/natmol);
572 pdba->atomname[n] = put_symtab(&symtab,
573 bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
574 if (bDiamond)
575 pdba->resname[n] = put_symtab(&symtab,"DIA");
576 else
577 pdba->resname[n] = put_symtab(&symtab,"SOL");
578 sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
579 pdba->atom[n].chain = ' ';
580 }
581
582 /* Generate the unit cell */
583 if (bDiamond)
584 unitcell_d(xx,box,odist);
585 else if (opt2bSet("-f",NFILE,fnm)) {
586 random_h_coords(natmol,natom/natmol,xx,box,bYaw,odist,hdist);
587 }
588 else
589 unitcell(xx,box,bYaw,odist,hdist);
590 if (debug) {
591 clear_mat(boxje);
592 boxje[XX][XX] = box[XX];
593 boxje[YY][YY] = box[YY];
594 boxje[ZZ][ZZ] = box[ZZ];
595 }
596 n=0;
597 for(i=0; (i<nx); i++) {
598 tmp[XX] = i*box[XX];
599 for(j=0; (j<ny); j++) {
600 tmp[YY] = j*box[YY];
601 for(k=0; (k<nz); k++) {
602 tmp[ZZ] = k*box[ZZ];
603 for(m=0; (m<natom); m++,n++) {
604 if ((!bOrdered && ((m % natmol) == 0)) || bOrdered)
605 rvec_add(xx[n % natom],tmp,xx[n]);
606 else
607 ;
608 }
609 }
610 }
611 }
612
613 clear_mat(boxje);
614 boxje[XX][XX] = box[XX]*nx;
615 boxje[YY][YY] = box[YY]*ny;
616 boxje[ZZ][ZZ] = box[ZZ]*nz;
617
618 printf("Crystal: %10.5f %10.5f %10.5f\n",
619 nx*box[XX],ny*box[YY],nz*box[ZZ]);
620
621 if (debug && !bDiamond) {
622 if (bSeries)
623 for(i=3; (i<=10*rcut); i++) {
624 fprintf(debug,"This is with rcut = %g\n",i*0.1);
625 virial(debug,bFull,nmax/natmol,xx,boxje,
626 0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
627 }
628 else
629 virial(debug,bFull,nmax/natmol,xx,boxje,
630 rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
631 }
632
633 if (bDiamond)
634 mk_diamond(pdba,xx,odist,&symtab,bPBC,boxje);
635
636 fn = ftp2fn(efSTO,NFILE,fnm);
637 if (fn2ftp(fn) == efPDB) {
638 fp = ffopen(fn,"w");
639 if (bDiamond)
640 fprintf(fp,"HEADER This is a *diamond*\n");
641 else
642 fprintf(fp,"HEADER A beautiful Ice Ih crystal\n");
643 fprintf(fp,"REMARK Generated by mkice with the following options:\n"
644 "REMARK nx = %d, ny = %d, nz = %d, odist = %g, hdist = %g\n",
645 nx,ny,nz,odist,hdist);
646 bromacs(quote,255);
647 write_pdbfile(fp,quote,pdba,xx,boxje,' ',-1);
648 ffclose(fp);
649 }
650 else {
651 bromacs(quote,255);
652 write_sto_conf(fn,quote,pdba,xx,NULL,boxje);
653 }
654
655 if (ftp2bSet(efTRN,NFILE,fnm))
656 write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
657
658 return 0;
659 }
660
The ultimate molecular dynamics simulation package