?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Partial commit of the project to remove all static variables.
[gromacs.git] / src / contrib / mkice.c
blobe6f3fe8a996b16bd2d2a8d66400c82cb5d72b6ef
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.1
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
16 *
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
23 *
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
26 *
27 * For more info, check our website at http://www.gromacs.org
28 *
29 * And Hey:
30 * Great Red Owns Many ACres of Sand
31 */
32
33 #include <stdio.h>
34 #include <math.h>
35 #include "typedefs.h"
36 #include "statutil.h"
37 #include "copyrite.h"
38 #include "fatal.h"
39 #include "pdbio.h"
40 #include "macros.h"
41 #include "smalloc.h"
42 #include "vec.h"
43 #include "pbc.h"
44 #include "physics.h"
45 #include "names.h"
46 #include "txtdump.h"
47 #include "trnio.h"
48 #include "symtab.h"
49 #include "atomprop.h"
50 #include "confio.h"
51
52 #define TET 109.47
53 #define DCONS 0.117265878
54
55 typedef struct {
56 int n,aa[4];
57 } t_bbb;
58
59 static char *watname[] = { "OW ", "HW1", "HW2", "DW", "SW" };
60 static char *diamname[] = { "C1", "H2", "H3", "H4", "H5", "H2", "H3", "H4", "H5" };
61 static real qspc[] = { -0.82, 0.41, 0.41 };
62 static real qyaw[] = { 1.24588, 0.62134, 0.62134, 0.0, -2.48856 };
63 static real spc_lj[3][6] = {
64 { 2.6171e-3, 2.6331e-6, 0, 0, 0, 0 },
65 { 0, 0, 0, 0, 0, 0 },
66 { 0, 0, 0, 0, 0, 0 }
67 };
68 #define CHH 2e-9
69 #define CHS 2e-9
70 #define COS 2e-6
71 static real yaw_lj[5][10] = {
72 { 0, 0, 0, 0, 0, 0, 0, 0, 0, COS },
73 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
74 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
75 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
76 { 0, COS, 0, CHS, 0, CHS, 0, 0, 2.6e-3, 0 }
77 };
78
79 void unitcell(rvec x[],rvec box,bool bYaw,real odist,real hdist)
80 {
81 #define cx 0.81649658
82 #define cy 0.47140452
83 #define cy2 0.94280904
84 #define cz 0.33333333
85
86 rvec xx[24] = {
87 { 0, 0, 0 }, /* O1 */
88 { 0, 0, 1 }, /* H relative to Oxygen */
89 { cx, cy, -cz },
90 { cx, cy, -cz }, /* O2 */
91 { 0, 0, -1 }, /* H relative to Oxygen */
92 { cx,-cy, +cz },
93 { cx, cy+cy2, 0 }, /* O3 */
94 { -cx, cy, -cz }, /* H relative to Oxygen */
95 { 0, -cy2, -cz },
96 { 0, 2*cy+cy2, -cz }, /* O4 */
97 {-cx,-cy, +cz }, /* H relative to Oxygen */
98 { 0 , cy2, +cz },
99 { 0, 0, 1 }, /* O5 */
100 {-cx, cy, +cz }, /* H relative to Oxygen */
101 { 0 , -cy2, +cz },
102 { cx, cy, 1+cz }, /* O6 */
103 { -cx, -cy, -cz }, /* H relative to Oxygen */
104 { 0, cy2, -cz },
105 { cx, cy+cy2, 1 }, /* O7 */
106 { 0, 0, -1 }, /* H relative to Oxygen */
107 { cx, cy, +cz },
108 { 0, 2*cy+cy2,1+cz }, /* O8 */
109 { 0, 0, 1 }, /* H relative to Oxygen */
110 { cx, -cy, -cz }
111 };
112 int i,iin,iout,j,m;
113 rvec tmp,t2,dip;
114
115 clear_rvec(dip);
116 for(i=0; (i<8); i++) {
117 iin = 3*i;
118 if (bYaw)
119 iout = 5*i;
120 else
121 iout = iin;
122 svmul(odist,xx[iin],x[iout]);
123 svmul(-0.82,x[iout],t2);
124 rvec_inc(dip,t2);
125 for(j=1; (j<=2); j++) {
126 svmul(hdist,xx[iin+j],tmp);
127 rvec_add(x[iout],tmp,x[iout+j]);
128 svmul(0.41,x[iout+j],t2);
129 rvec_inc(dip,t2);
130 }
131 if (bYaw)
132 for(m=0; (m<DIM); m++)
133 x[iout+3][m] = x[iout+4][m] =
134 (1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
135 }
136
137 box[XX] = 2*cx;
138 box[YY] = 2*(cy2+cy);
139 box[ZZ] = 2*(1+cz);
140 for(i=0; (i<DIM); i++)
141 box[i] *= odist;
142
143 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
144 printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
145 }
146
147 void unitcell_d(rvec x[],rvec box,real odist)
148 {
149 rvec cc[8] = {
150 { 0, 0, 0 }, /* C1 */
151 { cx, cy, -cz }, /* C2 */
152 { cx, cy+cy2, 0 }, /* C3 */
153 { 0, 2*cy+cy2, -cz }, /* C4 */
154 { 0, 0, 1 }, /* C5 */
155 { cx, cy, 1+cz }, /* C6 */
156 { cx, cy+cy2, 1 }, /* C7 */
157 { 0, 2*cy+cy2,1+cz }, /* C8 */
158 };
159 rvec hh[4] = {
160 { 0, 0, 1 }, /* H relative to C */
161 { cx, cy, -cz },
162 { cx, -cy, -cz },
163 {-cy2, 0, -cz }
164 };
165 int i,iin,iout,j,m;
166 rvec tmp,t2,dip;
167
168 clear_rvec(dip);
169 for(i=0; (i<8); i++) {
170 iin = i;
171 iout = i;
172 svmul(odist,cc[iin],x[iout]);
173 }
174 box[XX] = 2*cx;
175 box[YY] = 2*(cy2+cy);
176 box[ZZ] = 2*(1+cz);
177 for(i=0; (i<DIM); i++)
178 box[i] *= odist;
179
180 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
181 }
182
183 static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
184 bool bPBC,matrix box)
185 {
186 real od2 = odist*odist+1e-5;
187 t_bbb *bbb;
188 int i,j;
189 rvec dx;
190
191 if (bPBC)
192 init_pbc(box);
193 snew(bbb,natoms);
194 for(i=0; (i<natoms); i++) {
195 for(j=i+1; (j<natoms); j++) {
196 if (bPBC)
197 pbc_dx(x[i],x[j],dx);
198 else
199 rvec_sub(x[i],x[j],dx);
200 if (iprod(dx,dx) <= od2) {
201 bbb[i].aa[bbb[i].n++] = j;
202 bbb[j].aa[bbb[j].n++] = i;
203 }
204 }
205 }
206 if (debug)
207 #define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
208 for(i=0; (i<natoms); i++)
209 fprintf(debug,"bonds from %3d: %d %d %d %d\n",
210 i,PRB(1),PRB(2),PRB(3),PRB(4));
211 #undef PRB
212 return bbb;
213 }
214
215 static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
216 bool bPBC,matrix box)
217 {
218
219 int i,ib,j,k,l,m,nrm=0;
220 t_bbb *bbb;
221 bool *bRemove;
222 rvec dx;
223
224 do {
225 nrm = 0;
226 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
227
228 for(i=0; (i<a->nr); i++) {
229 if (bbb[i].n < 2) {
230 for(k=0; (k<bbb[i].n); k++) {
231 ib = bbb[i].aa[k];
232 for(j=0; (j<bbb[ib].n); j++)
233 if (bbb[ib].aa[j] == i)
234 break;
235 if (j == bbb[ib].n)
236 fatal_error(0,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
237 for( ; (j<bbb[ib].n-1); j++)
238 bbb[ib].aa[j] = bbb[ib].aa[j+1];
239 bbb[ib].n--;
240 nrm++;
241 }
242 bbb[i].n = 0;
243 }
244 }
245
246 for(i=j=0; (i<a->nr); i++) {
247 if (bbb[i].n >= 2) {
248 copy_rvec(x[i],x[j]);
249 j++;
250 }
251 }
252 fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
253 a->nr-j,a->nr);
254 a->nr = j;
255 sfree(bbb);
256 } while (nrm > 0);
257 /* Rename atoms */
258 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
259 for(i=0; (i<a->nr); i++) {
260 switch (bbb[i].n) {
261 case 4:
262 a->atomname[i] = put_symtab(symtab,"C");
263 break;
264 case 3:
265 a->atomname[i] = put_symtab(symtab,"CH1");
266 break;
267 case 2:
268 a->atomname[i] = put_symtab(symtab,"CH2");
269 break;
270 default:
271 fatal_error(0,"This atom (%d) has %d bonds only",i,bbb[i].n);
272 }
273 }
274 sfree(bbb);
275 }
276
277 static real calc_ener(real c6,real c12,rvec dx,tensor vir)
278 {
279 real r,e,f;
280 int m,n;
281
282 r = norm(dx);
283 e = c12*pow(r,-12.0) - c6*pow(r,-6.0);
284 f = 12*c12*pow(r,-14.0) - 6*c6*pow(r,-8.0);
285 for(m=0; (m<DIM); m++)
286 for(n=0; (n<DIM); n++)
287 vir[m][n] -= 0.5*f*dx[m]*dx[n];
288
289 return e;
290 }
291
292 void virial(FILE *fp,bool bFull,int nmol,rvec x[],matrix box,real rcut,
293 bool bYaw,real q[],bool bLJ)
294 {
295 int i,j,im,jm,natmol,ik,jk,m,ninter;
296 rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
297 real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
298
299 init_pbc(box);
300 fprintf(fp,"%3s - %3s: %6s %6s %6s %6s %8s %8s %8s\n",
301 "ai","aj","dx","dy","dz","|d|","virx","viry","virz");
302 clear_rvec(ftot);
303 clear_rvec(vir);
304 ninter = 0;
305 vctot = 0;
306 v12tot = 0;
307 v6tot = 0;
308 natmol = bYaw ? 5 : 3;
309 snew(force,nmol*natmol);
310
311 for(i=0; (i<nmol); i++) {
312 im = natmol*i;
313 /* Center of geometry */
314 clear_rvec(xcmi);
315 for(ik=0; (ik<natmol); ik++)
316 rvec_inc(xcmi,x[im+ik]);
317 for(m=0; (m<DIM); m++)
318 xcmi[m] /= natmol;
319
320 for(j=i+1; (j<nmol); j++) {
321 jm = natmol*j;
322 /* Center of geometry */
323 clear_rvec(xcmj);
324 for(jk=0; (jk<natmol); jk++)
325 rvec_inc(xcmj,x[jm+jk]);
326 for(m=0; (m<DIM); m++)
327 xcmj[m] /= natmol;
328
329 /* First check COM-COM distance */
330 pbc_dx(xcmi,xcmj,dx);
331 if (norm(dx) < rcut) {
332 ninter++;
333 /* Neirest neighbour molecules! */
334 clear_rvec(dvir);
335 for(ik=0; (ik<natmol); ik++) {
336 for(jk=0; (jk<natmol); jk++) {
337 pbc_dx(x[im+ik],x[jm+jk],dx);
338 dx2 = iprod(dx,dx);
339 dx1 = sqrt(dx2);
340 vcoul = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
341 vctot += vcoul;
342
343 if (bLJ) {
344 if (bYaw) {
345 c6 = yaw_lj[ik][2*jk];
346 c12 = yaw_lj[ik][2*jk+1];
347 }
348 else {
349 c6 = spc_lj[ik][2*jk];
350 c12 = spc_lj[ik][2*jk+1];
351 }
352 dx6 = dx2*dx2*dx2;
353 v6 = c6/dx6;
354 v12 = c12/(dx6*dx6);
355 v6tot -= v6;
356 v12tot+= v12;
357 fscal = (vcoul+12*v12-6*v6)/dx2;
358 }
359 else
360 fscal = vcoul/dx2;
361 for(m=0; (m<DIM); m++) {
362 f[m] = dx[m]*fscal;
363 dvir[m] -= 0.5*dx[m]*f[m];
364 }
365 rvec_inc(force[ik+im],f);
366 rvec_dec(force[jk+jm],f);
367 /*if (bFull)
368 fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
369 " %8.3f %8.3f %8.3f\n",
370 watname[ik],im+ik,watname[jk],jm+jk,
371 dx[XX],dx[YY],dx[ZZ],norm(dx),
372 dvir[XX],dvir[YY],dvir[ZZ]);*/
373 }
374 }
375 if (bFull)
376 fprintf(fp,"%3s%4d-%3s%4d: "
377 " %8.3f %8.3f %8.3f\n",
378 "SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
379 rvec_inc(vir,dvir);
380 }
381 }
382 }
383 fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
384 ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
385 fprintf(fp,"Vcoul: %10.4e V12: %10.4e V6: %10.4e Vtot: %10.4e (kJ/mol)\n",
386 vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
387 pr_rvec(fp,0,"vir ",vir,DIM);
388
389 for(m=0; (m<DIM); m++)
390 pres[m] = -2*PRESFAC/(det(box))*vir[m];
391 pr_rvec(fp,0,"pres",pres,DIM);
392 pr_rvecs(fp,0,"force",force,natmol*nmol);
393 sfree(force);
394 }
395
396 static real calc_mass(t_atoms *atoms,bool bGetMass)
397 {
398 real tmass;
399 int i;
400
401 tmass = 0;
402 for(i=0; (i<atoms->nr); i++) {
403 if (bGetMass)
404 atoms->atom[i].m = get_mass(*atoms->resname[atoms->atom[i].resnr],
405 *atoms->atomname[i]);
406 tmass += atoms->atom[i].m;
407 }
408
409 return tmass;
410 }
411
412 static real density(t_atoms *atoms,matrix box)
413 {
414 real tm;
415 real vol;
416
417 vol = det(box);
418 tm = calc_mass(atoms,TRUE);
419
420 return (tm*AMU)/(vol*NANO*NANO*NANO);
421 }
422
423 int main(int argc,char *argv[])
424 {
425 static char *desc[] = {
426 "mkice generates an ice crystal in the Ih crystal form which is the",
427 "most stable form. The rectangular unitcell contains eight molecules",
428 "and all oxygens are tetrahedrally coordinated."
429 };
430 static int nx=1,ny=1,nz=1;
431 static bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE;
432 static bool bOrdered=TRUE,bDiamond=FALSE,bPBC=TRUE;
433 static real rcut=0.3,odist=0.274,hdist=0.09572;
434 t_pargs pa[] = {
435 { "-nx", FALSE, etINT, {&nx}, "nx" },
436 { "-ny", FALSE, etINT, {&ny}, "ny" },
437 { "-nz", FALSE, etINT, {&nz}, "nz" },
438 { "-yaw", FALSE, etBOOL, {&bYaw},
439 "Generate dummies and shell positions" },
440 { "-lj", FALSE, etBOOL, {&bLJ},
441 "Use LJ as well as coulomb for virial calculation" },
442 { "-rcut", FALSE,etREAL, {&rcut},"Cut-off for virial calculations" },
443 { "-full", FALSE,etBOOL, {&bFull},"Full virial output" },
444 { "-odist", FALSE, etREAL, {&odist}, "Distance (nm) between oxygens" },
445 { "-hdist", FALSE, etREAL, {&hdist}, "Bondlength (nm) for OH bond" },
446 { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
447 { "-pbc", FALSE, etBOOL, {&bPBC}, "Make a pariodic diamond" },
448 { "-order", FALSE,etBOOL, {&bOrdered}, "Make a proton-ordered ice lattice" },
449 { "-series",FALSE, etBOOL, {&bSeries},
450 "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
451 };
452 t_filenm fnm[] = {
453 { efSTO, "-p", "ice", ffWRITE },
454 { efSTO, "-c", NULL, ffOPTRD },
455 { efTRN, "-o", "ice", ffOPTWR }
456 };
457 #define NFILE asize(fnm)
458
459 FILE *fp;
460 char *fn,quote[256];
461 int i,j,k,n,nmax,m,natom,natmol;
462 t_atoms *pdba;
463 t_atoms atoms;
464 t_symtab symtab;
465 rvec box,tmp,*xx;
466 matrix boxje;
467
468 CopyRight(stdout,argv[0]);
469 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),
470 desc,0,NULL);
471 if (debug) {
472 fprintf(debug,"nx = %3d, ny = %3d, nz = %3d\n",nx,ny,nz);
473 fprintf(debug,"YAW = %3s, LJ = %3s, rcut = %g\n",yesno_names[bYaw],
474 yesno_names[bLJ],rcut);
475 }
476
477 if (bYaw)
478 natmol = 5;
479 else if (bDiamond)
480 natmol = 1;
481 else
482 natmol = 3;
483 natom = natmol*8;
484 nmax = natom*nx*ny*nz;
485 snew(xx,nmax);
486 snew(pdba,1);
487 init_t_atoms(pdba,nmax,TRUE);
488 pdba->nr = nmax;
489 open_symtab(&symtab);
490 for(n=0; (n<nmax); n++) {
491 pdba->pdbinfo[n].type = epdbATOM;
492 pdba->pdbinfo[n].atomnr = n;
493 pdba->atom[n].resnr = n/natmol;
494 pdba->atomname[n] = put_symtab(&symtab,
495 bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
496 if (bDiamond)
497 pdba->resname[n] = put_symtab(&symtab,"DIA");
498 else
499 pdba->resname[n] = put_symtab(&symtab,"SOL");
500 sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
501 pdba->atom[n].chain = ' ';
502 }
503
504 /* Generate the unit cell */
505 if (bDiamond)
506 unitcell_d(xx,box,odist);
507 else
508 unitcell(xx,box,bYaw,odist,hdist);
509 if (debug) {
510 clear_mat(boxje);
511 boxje[XX][XX] = box[XX];
512 boxje[YY][YY] = box[YY];
513 boxje[ZZ][ZZ] = box[ZZ];
514 }
515 n=0;
516 for(i=0; (i<nx); i++) {
517 tmp[XX] = i*box[XX];
518 for(j=0; (j<ny); j++) {
519 tmp[YY] = j*box[YY];
520 for(k=0; (k<nz); k++) {
521 tmp[ZZ] = k*box[ZZ];
522 for(m=0; (m<natom); m++,n++) {
523 if ((!bOrdered && ((m % natmol) == 0)) || bOrdered)
524 rvec_add(xx[n % natom],tmp,xx[n]);
525 else
526 ;
527 }
528 }
529 }
530 }
531
532 clear_mat(boxje);
533 boxje[XX][XX] = box[XX]*nx;
534 boxje[YY][YY] = box[YY]*ny;
535 boxje[ZZ][ZZ] = box[ZZ]*nz;
536
537 printf("Crystal: %10.5f %10.5f %10.5f\n",
538 nx*box[XX],ny*box[YY],nz*box[ZZ]);
539
540 if (debug && !bDiamond) {
541 if (bSeries)
542 for(i=3; (i<=10*rcut); i++) {
543 fprintf(debug,"This is with rcut = %g\n",i*0.1);
544 virial(debug,bFull,nmax/natmol,xx,boxje,
545 0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
546 }
547 else
548 virial(debug,bFull,nmax/natmol,xx,boxje,
549 rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
550 }
551
552 if (bDiamond)
553 mk_diamond(pdba,xx,odist,&symtab,bPBC,boxje);
554
555 fn = ftp2fn(efSTO,NFILE,fnm);
556 if (fn2ftp(fn) == efPDB) {
557 fp = ffopen(fn,"w");
558 if (bDiamond)
559 fprintf(fp,"HEADER This is a *diamond*\n");
560 else
561 fprintf(fp,"HEADER A beautiful Ice Ih crystal\n");
562 fprintf(fp,"REMARK Generated by mkice with the following options:\n"
563 "REMARK nx = %d, ny = %d, nz = %d, odist = %g, hdist = %g\n",
564 nx,ny,nz,odist,hdist);
565 fprintf(fp,"REMARK Density of this crystal is %g (g/l)\n",
566 density(pdba,boxje));
567 write_pdbfile(fp,bromacs(quote,255),pdba,xx,boxje,' ',-1);
568 fclose(fp);
569 }
570 else {
571 write_sto_conf(fn,bromacs(quote,255),pdba,xx,NULL,boxje);
572 }
573
574 if (ftp2bSet(efTRN,NFILE,fnm))
575 write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
576
577 return 0;
578 }
579
The ultimate molecular dynamics simulation package