moldyn.c

   1 /*
   2  * moldyn.c - molecular dynamics library main file
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 #include "moldyn.h"
  19
  20 #include "math/math.h"
  21 #include "init/init.h"
  22 #include "random/random.h"
  23 #include "visual/visual.h"
  24 #include "list/list.h"
  25
  26 int moldyn_usage(char **argv) {
  27
  28         printf("\n%s usage:\n\n",argv[0]);
  29         printf("--- general options ---\n");
  30         printf("-E <steps> <file> (log total energy)\n");
  31         printf("-M <steps> <file> (log total momentum)\n");
  32         printf("-D <steps> <file> (dump total information)\n");
  33         printf("-S <steps> <filebase> (single save file)\n");
  34         printf("-V <steps> <filebase> (rasmol file)\n");
  35         printf("--- physics options ---\n");
  36         printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
  37         printf("-t <timestep tau> [s] (%.15f)\n",MOLDYN_TAU);
  38         printf("-C <cutoff radius> [m] (%.15f)\n",MOLDYN_CUTOFF);
  39         printf("-R <runs> (%d)\n",MOLDYN_RUNS);
  40         printf(" -- integration algo --\n");
  41         printf("  -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
  42         printf("     0: velocity verlet\n");
  43         printf(" -- potential --\n");
  44         printf("  -P <number> <param1 param2 ...>\n");
  45         printf("     0: harmonic oscillator\n");
  46         printf("        param1: spring constant\n");
  47         printf("        param2: equilibrium distance\n");
  48         printf("     1: lennard jones\n");
  49         printf("        param1: epsilon\n");
  50         printf("        param2: sigma\n");
  51         printf("\n");
  52
  53         return 0;
  54 }
  55
  56 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
  57
  58         int i;
  59         t_ho_params hop;
  60         t_lj_params ljp;
  61         double s,e;
  62
  63         memset(moldyn,0,sizeof(t_moldyn));
  64
  65         /* default values */
  66         moldyn->t=MOLDYN_TEMP;
  67         moldyn->tau=MOLDYN_TAU;
  68         moldyn->time_steps=MOLDYN_RUNS;
  69         moldyn->integrate=velocity_verlet;
  70         moldyn->potential_force_function=lennard_jones;
  71
  72         /* parse argv */
  73         for(i=1;i<argc;i++) {
  74                 if(argv[i][0]=='-') {
  75                         switch(argv[i][1]){
  76                                 case 'E':
  77                                         moldyn->ewrite=atoi(argv[++i]);
  78                                         strncpy(moldyn->efb,argv[++i],64);
  79                                         break;
  80                                 case 'M':
  81                                         moldyn->mwrite=atoi(argv[++i]);
  82                                         strncpy(moldyn->mfb,argv[++i],64);
  83                                         break;
  84                                 case 'D':
  85                                         moldyn->dwrite=atoi(argv[++i]);
  86                                         strncpy(moldyn->dfb,argv[++i],64);
  87                                         break;
  88                                 case 'S':
  89                                         moldyn->swrite=atoi(argv[++i]);
  90                                         strncpy(moldyn->sfb,argv[++i],64);
  91                                         break;
  92                                 case 'V':
  93                                         moldyn->vwrite=atoi(argv[++i]);
  94                                         strncpy(moldyn->vfb,argv[++i],64);
  95                                         break;
  96                                 case 'T':
  97                                         moldyn->t=atof(argv[++i]);
  98                                         break;
  99                                 case 't':
 100                                         moldyn->tau=atof(argv[++i]);
 101                                         break;
 102                                 case 'C':
 103                                         moldyn->cutoff=atof(argv[++i]);
 104                                         break;
 105                                 case 'R':
 106                                         moldyn->time_steps=atoi(argv[++i]);
 107                                         break;
 108                                 case 'I':
 109         /* integration algorithm */
 110         switch(atoi(argv[++i])) {
 111                 case MOLDYN_INTEGRATE_VERLET:
 112                         moldyn->integrate=velocity_verlet;
 113                         break;
 114                 default:
 115                         printf("unknown integration algo %s\n",argv[i]);
 116                         moldyn_usage(argv);
 117                         return -1;
 118         }
 119
 120                                 case 'P':
 121         /* potential + params */
 122         switch(atoi(argv[++i])) {
 123                 case MOLDYN_POTENTIAL_HO:
 124                         hop.spring_constant=atof(argv[++i]);
 125                         hop.equilibrium_distance=atof(argv[++i]);
 126                         moldyn->pot_params=malloc(sizeof(t_ho_params));
 127                         memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
 128                         moldyn->potential_force_function=harmonic_oscillator;
 129                         break;
 130                 case MOLDYN_POTENTIAL_LJ:
 131                         e=atof(argv[++i]);
 132                         s=atof(argv[++i]);
 133                         ljp.epsilon4=4*e;
 134                         ljp.sigma6=s*s*s*s*s*s;
 135                         ljp.sigma12=ljp.sigma6*ljp.sigma6;
 136                         moldyn->pot_params=malloc(sizeof(t_lj_params));
 137                         memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
 138                         moldyn->potential_force_function=lennard_jones;
 139                         break;
 140                 default:
 141                         printf("unknown potential %s\n",argv[i]);
 142                         moldyn_usage(argv);
 143                         return -1;
 144         }
 145
 146                                 default:
 147                                         printf("unknown option %s\n",argv[i]);
 148                                         moldyn_usage(argv);
 149                                         return -1;
 150                         }
 151                 } else {
 152                         moldyn_usage(argv);
 153                         return -1;
 154                 }
 155         }
 156
 157         return 0;
 158 }
 159
 160 int moldyn_log_init(t_moldyn *moldyn) {
 161
 162         moldyn->lvstat=0;
 163         t_visual *vis;
 164
 165         vis=&(moldyn->vis);
 166
 167         if(moldyn->ewrite) {
 168                 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
 169                 if(moldyn->efd<0) {
 170                         perror("[moldyn] efd open");
 171                         return moldyn->efd;
 172                 }
 173                 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
 174                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
 175         }
 176
 177         if(moldyn->mwrite) {
 178                 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
 179                 if(moldyn->mfd<0) {
 180                         perror("[moldyn] mfd open");
 181                         return moldyn->mfd;
 182                 }
 183                 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
 184                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
 185         }
 186
 187         if(moldyn->swrite)
 188                 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
 189
 190         if(moldyn->dwrite) {
 191                 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
 192                 if(moldyn->dfd<0) {
 193                         perror("[moldyn] dfd open");
 194                         return moldyn->dfd;
 195                 }
 196                 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
 197                 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
 198         }
 199
 200         if((moldyn->vwrite)&&(vis)) {
 201                 moldyn->visual=vis;
 202                 visual_init(vis,moldyn->vfb);
 203                 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
 204         }
 205
 206         moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
 207
 208         return 0;
 209 }
 210
 211 int moldyn_log_shutdown(t_moldyn *moldyn) {
 212
 213         if(moldyn->efd) close(moldyn->efd);
 214         if(moldyn->mfd) close(moldyn->efd);
 215         if(moldyn->dfd) close(moldyn->efd);
 216         if(moldyn->visual) visual_tini(moldyn->visual);
 217
 218         return 0;
 219 }
 220
 221 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 222
 223         int ret;
 224
 225         ret=moldyn_parse_argv(moldyn,argc,argv);
 226         if(ret<0) return ret;
 227
 228         ret=moldyn_log_init(moldyn);
 229         if(ret<0) return ret;
 230
 231         rand_init(&(moldyn->random),NULL,1);
 232         moldyn->random.status|=RAND_STAT_VERBOSE;
 233
 234         moldyn->status=0;
 235
 236         return 0;
 237 }
 238
 239 int moldyn_shutdown(t_moldyn *moldyn) {
 240
 241         moldyn_log_shutdown(moldyn);
 242         rand_close(&(moldyn->random));
 243         free(moldyn->atom);
 244
 245         return 0;
 246 }
 247
 248 int create_lattice(unsigned char type,int element,double mass,double lc,
 249                    int a,int b,int c,t_atom **atom) {
 250
 251         int count;
 252         int ret;
 253         t_3dvec origin;
 254
 255         count=a*b*c;
 256
 257         if(type==FCC) count*=4;
 258         if(type==DIAMOND) count*=8;
 259
 260         *atom=malloc(count*sizeof(t_atom));
 261         if(*atom==NULL) {
 262                 perror("malloc (atoms)");
 263                 return -1;
 264         }
 265
 266         v3_zero(&origin);
 267
 268         switch(type) {
 269                 case FCC:
 270                         ret=fcc_init(a,b,c,lc,*atom,&origin);
 271                         break;
 272                 case DIAMOND:
 273                         ret=diamond_init(a,b,c,lc,*atom,&origin);
 274                         break;
 275                 default:
 276                         printf("unknown lattice type (%02x)\n",type);
 277                         return -1;
 278         }
 279
 280         /* debug */
 281         if(ret!=count) {
 282                 printf("ok, there is something wrong ...\n");
 283                 printf("calculated -> %d atoms\n",count);
 284                 printf("created -> %d atoms\n",ret);
 285                 return -1;
 286         }
 287
 288         while(count) {
 289                 (*atom)[count-1].element=element;
 290                 (*atom)[count-1].mass=mass;
 291                 count-=1;
 292         }
 293
 294         return ret;
 295 }
 296
 297 int destroy_lattice(t_atom *atom) {
 298
 299         if(atom) free(atom);
 300
 301         return 0;
 302 }
 303
 304 int thermal_init(t_moldyn *moldyn) {
 305
 306         /*
 307          * - gaussian distribution of velocities
 308          * - zero total momentum
 309          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 310          */
 311
 312         int i;
 313         double v,sigma;
 314         t_3dvec p_total,delta;
 315         t_atom *atom;
 316         t_random *random;
 317
 318         atom=moldyn->atom;
 319         random=&(moldyn->random);
 320
 321         /* gaussian distribution of velocities */
 322         v3_zero(&p_total);
 323         for(i=0;i<moldyn->count;i++) {
 324                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
 325                 /* x direction */
 326                 v=sigma*rand_get_gauss(random);
 327                 atom[i].v.x=v;
 328                 p_total.x+=atom[i].mass*v;
 329                 /* y direction */
 330                 v=sigma*rand_get_gauss(random);
 331                 atom[i].v.y=v;
 332                 p_total.y+=atom[i].mass*v;
 333                 /* z direction */
 334                 v=sigma*rand_get_gauss(random);
 335                 atom[i].v.z=v;
 336                 p_total.z+=atom[i].mass*v;
 337         }
 338
 339         /* zero total momentum */
 340         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 341         for(i=0;i<moldyn->count;i++) {
 342                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 343                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 344         }
 345
 346         /* velocity scaling */
 347         scale_velocity(moldyn);
 348
 349         return 0;
 350 }
 351
 352 int scale_velocity(t_moldyn *moldyn) {
 353
 354         int i;
 355         double e,c;
 356         t_atom *atom;
 357
 358         atom=moldyn->atom;
 359
 360         /*
 361          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 362          */
 363         e=0.0;
 364         for(i=0;i<moldyn->count;i++)
 365                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 366         c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
 367         for(i=0;i<moldyn->count;i++)
 368                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 369
 370         return 0;
 371 }
 372
 373 double get_e_kin(t_atom *atom,int count) {
 374
 375         int i;
 376         double e;
 377
 378         e=0.0;
 379
 380         for(i=0;i<count;i++) {
 381                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 382         }
 383
 384         return e;
 385 }
 386
 387 double get_e_pot(t_moldyn *moldyn) {
 388
 389         return moldyn->energy;
 390 }
 391
 392 double get_total_energy(t_moldyn *moldyn) {
 393
 394         double e;
 395
 396         e=get_e_kin(moldyn->atom,moldyn->count);
 397         e+=get_e_pot(moldyn);
 398
 399         return e;
 400 }
 401
 402 t_3dvec get_total_p(t_atom *atom, int count) {
 403
 404         t_3dvec p,p_total;
 405         int i;
 406
 407         v3_zero(&p_total);
 408         for(i=0;i<count;i++) {
 409                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 410                 v3_add(&p_total,&p_total,&p);
 411         }
 412
 413         return p_total;
 414 }
 415
 416 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
 417
 418         double tau;
 419
 420         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
 421         tau*=1.0E-9;
 422         if(tau<moldyn->tau)
 423                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
 424                        moldyn->tau,tau);
 425
 426         return tau;
 427 }
 428
 429 /*
 430  * numerical tricks
 431  */
 432
 433 /* linked list / cell method */
 434
 435 int link_cell_init(t_moldyn *moldyn) {
 436
 437         t_linkcell *lc;
 438         int i;
 439
 440         lc=&(moldyn->lc);
 441
 442         /* list log fd */
 443         lc->listfd=open("/dev/null",O_WRONLY);
 444
 445         /* partitioning the md cell */
 446         lc->nx=moldyn->dim.x/moldyn->cutoff;
 447         lc->x=moldyn->dim.x/lc->nx;
 448         lc->ny=moldyn->dim.y/moldyn->cutoff;
 449         lc->y=moldyn->dim.y/lc->ny;
 450         lc->nz=moldyn->dim.z/moldyn->cutoff;
 451         lc->z=moldyn->dim.z/lc->nz;
 452
 453         lc->cells=lc->nx*lc->ny*lc->nz;
 454         lc->subcell=malloc(lc->cells*sizeof(t_list));
 455
 456         printf("initializing linked cells (%d)\n",lc->cells);
 457
 458         for(i=0;i<lc->cells;i++)
 459                 //list_init(&(lc->subcell[i]),1);
 460                 list_init(&(lc->subcell[i]));
 461
 462         link_cell_update(moldyn);
 463
 464         return 0;
 465 }
 466
 467 int link_cell_update(t_moldyn *moldyn) {
 468
 469         int count,i,j,k;
 470         int nx,ny,nz;
 471         t_atom *atom;
 472         t_linkcell *lc;
 473
 474         atom=moldyn->atom;
 475         lc=&(moldyn->lc);
 476
 477         nx=lc->nx;
 478         ny=lc->ny;
 479         nz=lc->nz;
 480
 481         for(i=0;i<lc->cells;i++)
 482                 list_destroy(&(moldyn->lc.subcell[i]));
 483
 484         for(count=0;count<moldyn->count;count++) {
 485                 i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x;
 486                 j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y;
 487                 k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z;
 488                 list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
 489                                        &(atom[count]));
 490         }
 491
 492         return 0;
 493 }
 494
 495 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 496
 497         t_linkcell *lc;
 498         int a;
 499         int count1,count2;
 500         int ci,cj,ck;
 501         int nx,ny,nz;
 502         int x,y,z;
 503         unsigned char bx,by,bz;
 504
 505         lc=&(moldyn->lc);
 506         nx=lc->nx;
 507         ny=lc->ny;
 508         nz=lc->nz;
 509         count1=1;
 510         count2=27;
 511         a=nx*ny;
 512
 513
 514         cell[0]=lc->subcell[i+j*nx+k*a];
 515         for(ci=-1;ci<=1;ci++) {
 516                 bx=0;
 517                 x=i+ci;
 518                 if((x<0)||(x>=nx)) {
 519                         x=(x+nx)%nx;
 520                         bx=1;
 521                 }
 522                 for(cj=-1;cj<=1;cj++) {
 523                         by=0;
 524                         y=j+cj;
 525                         if((y<0)||(y>=ny)) {
 526                                 y=(y+ny)%ny;
 527                                 by=1;
 528                         }
 529                         for(ck=-1;ck<=1;ck++) {
 530                                 bz=0;
 531                                 z=k+ck;
 532                                 if((z<0)||(z>=nz)) {
 533                                         z=(z+nz)%nz;
 534                                         bz=1;
 535                                 }
 536                                 if(!(ci|cj|ck)) continue;
 537                                 if(bx|by|bz) {
 538                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
 539                                 }
 540                                 else {
 541                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
 542                                 }
 543                         }
 544                 }
 545         }
 546
 547         return count2;
 548 }
 549
 550 int link_cell_shutdown(t_moldyn *moldyn) {
 551
 552         int i;
 553         t_linkcell *lc;
 554
 555         lc=&(moldyn->lc);
 556
 557         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
 558                 list_shutdown(&(moldyn->lc.subcell[i]));
 559
 560         if(lc->listfd) close(lc->listfd);
 561
 562         return 0;
 563 }
 564
 565 /*
 566  *
 567  * 'integration of newtons equation' - algorithms
 568  *
 569  */
 570
 571 /* start the integration */
 572
 573 int moldyn_integrate(t_moldyn *moldyn) {
 574
 575         int i;
 576         unsigned int e,m,s,d,v;
 577         t_3dvec p;
 578
 579         int fd;
 580         char fb[128];
 581
 582         /* initialize linked cell method */
 583         link_cell_init(moldyn);
 584
 585         /* logging & visualization */
 586         e=moldyn->ewrite;
 587         m=moldyn->mwrite;
 588         s=moldyn->swrite;
 589         d=moldyn->dwrite;
 590         v=moldyn->vwrite;
 591
 592         if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
 593                 printf("[moldyn] warning, lv system not initialized\n");
 594                 return -1;
 595         }
 596
 597         /* sqaure of some variables */
 598         moldyn->tau_square=moldyn->tau*moldyn->tau;
 599         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 600
 601         /* calculate initial forces */
 602         moldyn->potential_force_function(moldyn);
 603
 604         for(i=0;i<moldyn->time_steps;i++) {
 605
 606                 /* integration step */
 607                 moldyn->integrate(moldyn);
 608
 609                 /* check for log & visualization */
 610                 if(e) {
 611                         if(!(i%e))
 612                                 dprintf(moldyn->efd,
 613                                         "%.15f %.45f\n",i*moldyn->tau,
 614                                         get_total_energy(moldyn));
 615                 }
 616                 if(m) {
 617                         if(!(i%m)) {
 618                                 p=get_total_p(moldyn->atom,moldyn->count);
 619                                 dprintf(moldyn->mfd,
 620                                         "%.15f %.45f\n",i*moldyn->tau,
 621                                         v3_norm(&p));
 622                         }
 623                 }
 624                 if(s) {
 625                         if(!(i%s)) {
 626                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
 627                                          moldyn->t,i*moldyn->tau);
 628                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
 629                                 if(fd<0) perror("[moldyn] save fd open");
 630                                 else {
 631                                         write(fd,moldyn,sizeof(t_moldyn));
 632                                         write(fd,moldyn->atom,
 633                                               moldyn->count*sizeof(t_atom));
 634                                 }
 635                         }
 636                 }
 637                 if(d) {
 638                         if(!(i%d))
 639                                 write(moldyn->dfd,moldyn->atom,
 640                                       moldyn->count*sizeof(t_atom));
 641
 642                 }
 643                 if(v) {
 644                         if(!(i%v)) {
 645                                 visual_atoms(moldyn->visual,i*moldyn->tau,
 646                                              moldyn->atom,moldyn->count);
 647                                 printf("\rsteps: %d",i);
 648                                 fflush(stdout);
 649                         }
 650                 }
 651         }
 652
 653         return 0;
 654 }
 655
 656 /* velocity verlet */
 657
 658 int velocity_verlet(t_moldyn *moldyn) {
 659
 660         int i,count;
 661         double tau,tau_square;
 662         t_3dvec delta;
 663         t_atom *atom;
 664
 665         atom=moldyn->atom;
 666         count=moldyn->count;
 667         tau=moldyn->tau;
 668         tau_square=moldyn->tau_square;
 669
 670         for(i=0;i<count;i++) {
 671                 /* new positions */
 672                 v3_scale(&delta,&(atom[i].v),tau);
 673                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 674                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 675                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 676                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
 677
 678                 /* velocities */
 679                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 680                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 681         }
 682
 683         /* neighbour list update */
 684 printf("list update ...\n");
 685         link_cell_update(moldyn);
 686 printf("done\n");
 687
 688         /* forces depending on chosen potential */
 689 printf("calc potential/force ...\n");
 690         moldyn->potential_force_function(moldyn);
 691 printf("done\n");
 692
 693         for(i=0;i<count;i++) {
 694                 /* again velocities */
 695                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 696                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 697         }
 698
 699         return 0;
 700 }
 701
 702
 703 /*
 704  *
 705  * potentials & corresponding forces
 706  *
 707  */
 708
 709 /* harmonic oscillator potential and force */
 710
 711 int harmonic_oscillator(t_moldyn *moldyn) {
 712
 713         t_ho_params *params;
 714         t_atom *atom,*btom;
 715         t_linkcell *lc;
 716         t_list *this,neighbour[27];
 717         int i,j,c;
 718         int count;
 719         t_3dvec force,distance;
 720         double d,u;
 721         double sc,equi_dist;
 722         int ni,nj,nk;
 723
 724         params=moldyn->pot_params;
 725         atom=moldyn->atom;
 726         lc=&(moldyn->lc);
 727         sc=params->spring_constant;
 728         equi_dist=params->equilibrium_distance;
 729         count=moldyn->count;
 730
 731         /* reset energy counter */
 732         u=0.0;
 733
 734         for(i=0;i<count;i++) {
 735                 /* reset force */
 736                 v3_zero(&(atom[i].f));
 737
 738                 /* determine cell + neighbours */
 739                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
 740                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
 741                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
 742                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 743
 744                 /*
 745                  * processing cell of atom i
 746                  * => no need to check for empty list (1 element at minimum)
 747                  */
 748                 this=&(neighbour[0]);
 749                 list_reset(this);
 750                 do {
 751                         btom=this->current->data;
 752                         if(btom==&(atom[i]))
 753                                 continue;
 754                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 755                         d=v3_norm(&distance);
 756                         if(d<=moldyn->cutoff) {
 757                                 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 758                                 v3_scale(&force,&distance,
 759                                          -sc*(1.0-(equi_dist/d)));
 760                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 761                         }
 762                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 763
 764                 /*
 765                  * direct neighbour cells
 766                  * => no boundary condition check necessary
 767                  */
 768                 for(j=1;j<c;j++) {
 769                         this=&(neighbour[j]);
 770                         list_reset(this); /* there might not be a single atom */
 771                         if(this->start!=NULL) {
 772
 773                         do {
 774                                 btom=this->current->data;
 775                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 776                                 d=v3_norm(&distance);
 777                                 if(d<=moldyn->cutoff) {
 778                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 779                                         v3_scale(&force,&distance,
 780                                                  -sc*(1.0-(equi_dist/d)));
 781                                         v3_add(&(atom[i].f),&(atom[i].f),
 782                                                &force);
 783                                 }
 784                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 785
 786                         }
 787                 }
 788
 789                 /*
 790                  * indirect neighbour cells
 791                  * => check boundary conditions
 792                  */
 793                 for(j=c;j<27;j++) {
 794                         this=&(neighbour[j]);
 795                         list_reset(this); /* check boundary conditions */
 796                         if(this->start!=NULL) {
 797
 798                         do {
 799                                 btom=this->current->data;
 800                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 801                                 v3_per_bound(&distance,&(moldyn->dim));
 802                                 d=v3_norm(&distance);
 803                                 if(d<=moldyn->cutoff) {
 804                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 805                                         v3_scale(&force,&distance,
 806                                                  -sc*(1.0-(equi_dist/d)));
 807                                         v3_add(&(atom[i].f),&(atom[i].f),
 808                                                &force);
 809                                 }
 810                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 811
 812                         }
 813                 }
 814         }
 815
 816         moldyn->energy=0.5*u;
 817
 818         return 0;
 819 }
 820
 821 /* lennard jones potential & force for one sort of atoms */
 822
 823 int lennard_jones(t_moldyn *moldyn) {
 824
 825         t_lj_params *params;
 826         t_atom *atom,*btom;
 827         t_linkcell *lc;
 828         t_list *this,neighbour[27];
 829         int i,j,c;
 830         int count;
 831         t_3dvec force,distance;
 832         double d,h1,h2,u;
 833         double eps,sig6,sig12;
 834         int ni,nj,nk;
 835
 836         params=moldyn->pot_params;
 837         atom=moldyn->atom;
 838         lc=&(moldyn->lc);
 839         count=moldyn->count;
 840         eps=params->epsilon4;
 841         sig6=params->sigma6;
 842         sig12=params->sigma12;
 843
 844         /* reset energy counter */
 845         u=0.0;
 846
 847         for(i=0;i<count;i++) {
 848                 /* reset force */
 849                 v3_zero(&(atom[i].f));
 850
 851                 /* determine cell + neighbours */
 852                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
 853                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
 854                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
 855                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 856
 857                 /* processing cell of atom i */
 858                 this=&(neighbour[0]);
 859                 list_reset(this); /* list has 1 element at minimum */
 860                 do {
 861                         btom=this->current->data;
 862                         if(btom==&(atom[i]))
 863                                 continue;
 864                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 865                         d=v3_absolute_square(&distance);        /* 1/r^2 */
 866                         if(d<=moldyn->cutoff_square) {
 867                                 d=1.0/d;        /* 1/r^2 */
 868                                 h2=d*d;         /* 1/r^4 */
 869                                 h2*=d;          /* 1/r^6 */
 870                                 h1=h2*h2;       /* 1/r^12 */
 871                                 u+=eps*(sig12*h1-sig6*h2);
 872                                 h2*=d;          /* 1/r^8 */
 873                                 h1*=d;          /* 1/r^14 */
 874                                 h2*=6*sig6;
 875                                 h1*=12*sig12;
 876                                 d=+h1-h2;
 877                                 d*=eps;
 878                                 v3_scale(&force,&distance,d);
 879                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 880                         }
 881                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 882
 883                 /* neighbours not doing boundary condition overflow */
 884                 for(j=1;j<c;j++) {
 885                         this=&(neighbour[j]);
 886                         list_reset(this); /* there might not be a single atom */
 887                         if(this->start!=NULL) {
 888
 889                         do {
 890                                 btom=this->current->data;
 891                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 892                                 d=v3_absolute_square(&distance); /* r^2 */
 893                                 if(d<=moldyn->cutoff_square) {
 894                                         d=1.0/d;        /* 1/r^2 */
 895                                         h2=d*d;         /* 1/r^4 */
 896                                         h2*=d;          /* 1/r^6 */
 897                                         h1=h2*h2;       /* 1/r^12 */
 898                                         u+=eps*(sig12*h1-sig6*h2);
 899                                         h2*=d;          /* 1/r^8 */
 900                                         h1*=d;          /* 1/r^14 */
 901                                         h2*=6*sig6;
 902                                         h1*=12*sig12;
 903                                         d=+h1-h2;
 904                                         d*=eps;
 905                                         v3_scale(&force,&distance,d);
 906                                         v3_add(&(atom[i].f),&(atom[i].f),
 907                                                &force);
 908                                 }
 909                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 910
 911                         }
 912                 }
 913
 914                 /* neighbours due to boundary conditions */
 915                 for(j=c;j<27;j++) {
 916                         this=&(neighbour[j]);
 917                         list_reset(this); /* check boundary conditions */
 918                         if(this->start!=NULL) {
 919
 920                         do {
 921                                 btom=this->current->data;
 922                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 923                                 v3_per_bound(&distance,&(moldyn->dim));
 924                                 d=v3_absolute_square(&distance); /* r^2 */
 925                                 if(d<=moldyn->cutoff_square) {
 926                                         d=1.0/d;        /* 1/r^2 */
 927                                         h2=d*d;         /* 1/r^4 */
 928                                         h2*=d;          /* 1/r^6 */
 929                                         h1=h2*h2;       /* 1/r^12 */
 930                                         u+=eps*(sig12*h1-sig6*h2);
 931                                         h2*=d;          /* 1/r^8 */
 932                                         h1*=d;          /* 1/r^14 */
 933                                         h2*=6*sig6;
 934                                         h1*=12*sig12;
 935                                         d=+h1-h2;
 936                                         d*=eps;
 937                                         v3_scale(&force,&distance,d);
 938                                         v3_add(&(atom[i].f),&(atom[i].f),
 939                                                &force);
 940                                 }
 941                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 942
 943                         }
 944                 }
 945         }
 946
 947         moldyn->energy=0.5*u;
 948
 949         return 0;
 950 }
 951
 952 /* tersoff potential & force for 2 sorts of atoms */
 953
 954 int tersoff(t_moldyn *moldyn) {
 955
 956         t_tersoff_params *params;
 957         t_atom *atom,*btom,*ktom;
 958         t_linkcell *lc;
 959         t_list *this,*thisk,neighbour[27],neighbourk[27];
 960         int i,j,k,c,ck;
 961         int count;
 962         double u;
 963         int ni,nj,nk;
 964         int ki,kj,kk;
 965
 966
 967         params=moldyn->pot_params;
 968         atom=moldyn->atom;
 969         lc=&(moldyn->lc);
 970         count=moldyn->count;
 971
 972         /* reset enrgy counter */
 973         u=0.0;
 974
 975         for(i=0;i<count;i++) {
 976                 /* reset force */
 977                 v3_zero(&(atom[i].f));
 978
 979                 /* determin cell neighbours */
 980                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
 981                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
 982                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
 983                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 984
 985                 /*
 986                  * processing cell of atom i
 987                  * => no need to check for empty list (1 element at minimum)
 988                  */
 989                 this=&(neighbour[0]);
 990                 list_reset(this);
 991                 do {
 992                         btom=this->current->data;
 993                         if(btom==&(atom[i]))
 994                                 continue;
 995
 996                         /* 2 body stuff */
 997
 998                         /* determine cell neighbours of btom */
 999                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1000                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1001                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1002                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1003                                                      neighbourk);
1004
1005                         /* cell of btom */
1006                         thisk=&(neighbourk[0]);
1007                         list_reset(thisk);
1008                         do {
1009                                 ktom=thisk->current->data;
1010                                 if(ktom==btom)
1011                                         continue;
1012                                 if(ktom==&(atom[i]))
1013                                         continue;
1014
1015                                 /* 3 body stuff (1) */
1016
1017                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1018
1019                         /* direct neighbours of btom cell */
1020                         for(k=1;k<ck;k++) {
1021                                 thisk=&(neighbourk[k]);
1022                                 list_reset(thisk);
1023                                 if(thisk->start!=NULL) {
1024
1025                                 do {
1026                                         ktom=thisk->current->data;
1027                                         if(ktom==&(atom[i]))
1028                                                 continue;
1029
1030                                 /* 3 body stuff (2) */
1031
1032                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1033
1034                                 }
1035                         }
1036
1037                         /* indirect neighbours of btom cell */
1038                         for(k=ck;k<27;k++) {
1039                                 thisk=&(neighbourk[k]);
1040                                 list_reset(thisk);
1041                                 if(thisk->start!=NULL) {
1042
1043                                 do {
1044                                         ktom=thisk->current->data;
1045                                         if(ktom==&(atom[i]))
1046                                                 continue;
1047
1048                                 /* 3 body stuff */
1049
1050                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1051
1052                                 }
1053                         }
1054
1055
1056                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1057
1058                 /*
1059                  * direct neighbour cells of atom i
1060                  */
1061                 for(j=1;j<c;j++) {
1062                         this=&(neighbour[j]);
1063                         list_reset(this);
1064                         if(this->start!=NULL) {
1065
1066                         do {
1067                                 btom=this->current->data;
1068
1069                                 /* 2 body stuff */
1070
1071
1072                         /* determine cell neighbours of btom */
1073                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1074                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1075                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1076                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1077                                                      neighbourk);
1078
1079                         /* cell of btom */
1080                         thisk=&(neighbourk[0]);
1081                         list_reset(thisk);
1082                         do {
1083                                 ktom=thisk->current->data;
1084                                 if(ktom==btom)
1085                                         continue;
1086                                 if(ktom==&(atom[i]))
1087                                         continue;
1088
1089                                 /* 3 body stuff (1) */
1090
1091                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1092
1093                         /* direct neighbours of btom cell */
1094                         for(k=1;k<ck;k++) {
1095                                 thisk=&(neighbourk[k]);
1096                                 list_reset(thisk);
1097                                 if(thisk->start!=NULL) {
1098
1099                                 do {
1100                                         ktom=thisk->current->data;
1101                                         if(ktom==&(atom[i]))
1102                                                 continue;
1103
1104                                 /* 3 body stuff (2) */
1105
1106                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1107
1108                                 }
1109                         }
1110
1111                         /* indirect neighbours of btom cell */
1112                         for(k=ck;k<27;k++) {
1113                                 thisk=&(neighbourk[k]);
1114                                 list_reset(thisk);
1115                                 if(thisk->start!=NULL) {
1116
1117                                 do {
1118                                         ktom=thisk->current->data;
1119                                         if(ktom==&(atom[i]))
1120                                                 continue;
1121
1122                                 /* 3 body stuff (3) */
1123
1124                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1125
1126                                 }
1127                         }
1128
1129
1130                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1131
1132                         }
1133                 }
1134
1135                 /*
1136                  * indirect neighbour cells of atom i
1137                  */
1138                 for(j=c;j<27;j++) {
1139                         this=&(neighbour[j]);
1140                         list_reset(this);
1141                         if(this->start!=NULL) {
1142
1143                         do {
1144                                 btom=this->current->data;
1145
1146                                 /* 2 body stuff */
1147
1148
1149                         /* determine cell neighbours of btom */
1150                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1151                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1152                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1153                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1154                                                      neighbourk);
1155
1156                         /* cell of btom */
1157                         thisk=&(neighbourk[0]);
1158                         list_reset(thisk);
1159                         do {
1160                                 ktom=thisk->current->data;
1161                                 if(ktom==btom)
1162                                         continue;
1163                                 if(ktom==&(atom[i]))
1164                                         continue;
1165
1166                                 /* 3 body stuff (1) */
1167
1168                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1169
1170                         /* direct neighbours of btom cell */
1171                         for(k=1;k<ck;k++) {
1172                                 thisk=&(neighbourk[k]);
1173                                 list_reset(thisk);
1174                                 if(thisk->start!=NULL) {
1175
1176                                 do {
1177                                         ktom=thisk->current->data;
1178                                         if(ktom==&(atom[i]))
1179                                                 continue;
1180
1181                                 /* 3 body stuff (2) */
1182
1183                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1184
1185                                 }
1186                         }
1187
1188                         /* indirect neighbours of btom cell */
1189                         for(k=ck;k<27;k++) {
1190                                 thisk=&(neighbourk[k]);
1191                                 list_reset(thisk);
1192                                 if(thisk->start!=NULL) {
1193
1194                                 do {
1195                                         ktom=thisk->current->data;
1196                                         if(ktom==&(atom[i]))
1197                                                 continue;
1198
1199                                 /* 3 body stuff (3) */
1200
1201                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1202
1203                                 }
1204                         }
1205
1206
1207                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1208
1209                         }
1210                 }
1211
1212         }
1213
1214         moldyn->energy=0.5*u;
1215
1216         return 0;
1217 }
1218