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]),lc->listfd);
 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         link_cell_update(moldyn);
 685
 686         /* forces depending on chosen potential */
 687         moldyn->potential_force_function(moldyn);
 688
 689         for(i=0;i<count;i++) {
 690                 /* again velocities */
 691                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 692                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 693         }
 694
 695         return 0;
 696 }
 697
 698
 699 /*
 700  *
 701  * potentials & corresponding forces
 702  *
 703  */
 704
 705 /* harmonic oscillator potential and force */
 706
 707 int harmonic_oscillator(t_moldyn *moldyn) {
 708
 709         t_ho_params *params;
 710         t_atom *atom,*btom;
 711         t_linkcell *lc;
 712         t_list *this,neighbour[27];
 713         int i,j,c;
 714         int count;
 715         t_3dvec force,distance;
 716         double d,u;
 717         double sc,equi_dist;
 718         int ni,nj,nk;
 719
 720         params=moldyn->pot_params;
 721         atom=moldyn->atom;
 722         lc=&(moldyn->lc);
 723         sc=params->spring_constant;
 724         equi_dist=params->equilibrium_distance;
 725         count=moldyn->count;
 726
 727         /* reset energy counter */
 728         u=0.0;
 729
 730         for(i=0;i<count;i++) {
 731                 /* reset force */
 732                 v3_zero(&(atom[i].f));
 733
 734                 /* determine cell + neighbours */
 735                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
 736                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
 737                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
 738                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 739
 740                 /*
 741                  * processing cell of atom i
 742                  * => no need to check for empty list (1 element at minimum)
 743                  */
 744                 this=&(neighbour[0]);
 745                 list_reset(this);
 746                 do {
 747                         btom=this->current->data;
 748                         if(btom==&(atom[i]))
 749                                 continue;
 750                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 751                         d=v3_norm(&distance);
 752                         if(d<=moldyn->cutoff) {
 753                                 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 754                                 v3_scale(&force,&distance,
 755                                          -sc*(1.0-(equi_dist/d)));
 756                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 757                         }
 758                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 759
 760                 /*
 761                  * direct neighbour cells
 762                  * => no boundary condition check necessary
 763                  */
 764                 for(j=1;j<c;j++) {
 765                         this=&(neighbour[j]);
 766                         list_reset(this); /* there might not be a single atom */
 767                         if(this->start!=NULL) {
 768
 769                         do {
 770                                 btom=this->current->data;
 771                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 772                                 d=v3_norm(&distance);
 773                                 if(d<=moldyn->cutoff) {
 774                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 775                                         v3_scale(&force,&distance,
 776                                                  -sc*(1.0-(equi_dist/d)));
 777                                         v3_add(&(atom[i].f),&(atom[i].f),
 778                                                &force);
 779                                 }
 780                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 781
 782                         }
 783                 }
 784
 785                 /*
 786                  * indirect neighbour cells
 787                  * => check boundary conditions
 788                  */
 789                 for(j=c;j<27;j++) {
 790                         this=&(neighbour[j]);
 791                         list_reset(this); /* check boundary conditions */
 792                         if(this->start!=NULL) {
 793
 794                         do {
 795                                 btom=this->current->data;
 796                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 797                                 v3_per_bound(&distance,&(moldyn->dim));
 798                                 d=v3_norm(&distance);
 799                                 if(d<=moldyn->cutoff) {
 800                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 801                                         v3_scale(&force,&distance,
 802                                                  -sc*(1.0-(equi_dist/d)));
 803                                         v3_add(&(atom[i].f),&(atom[i].f),
 804                                                &force);
 805                                 }
 806                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 807
 808                         }
 809                 }
 810         }
 811
 812         moldyn->energy=0.5*u;
 813
 814         return 0;
 815 }
 816
 817 /* lennard jones potential & force for one sort of atoms */
 818
 819 int lennard_jones(t_moldyn *moldyn) {
 820
 821         t_lj_params *params;
 822         t_atom *atom,*btom;
 823         t_linkcell *lc;
 824         t_list *this,neighbour[27];
 825         int i,j,c;
 826         int count;
 827         t_3dvec force,distance;
 828         double d,h1,h2,u;
 829         double eps,sig6,sig12;
 830         int ni,nj,nk;
 831
 832         params=moldyn->pot_params;
 833         atom=moldyn->atom;
 834         lc=&(moldyn->lc);
 835         count=moldyn->count;
 836         eps=params->epsilon4;
 837         sig6=params->sigma6;
 838         sig12=params->sigma12;
 839
 840         /* reset energy counter */
 841         u=0.0;
 842
 843         for(i=0;i<count;i++) {
 844                 /* reset force */
 845                 v3_zero(&(atom[i].f));
 846
 847                 /* determine cell + neighbours */
 848                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
 849                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
 850                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
 851                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 852
 853                 /* processing cell of atom i */
 854                 this=&(neighbour[0]);
 855                 list_reset(this); /* list has 1 element at minimum */
 856                 do {
 857                         btom=this->current->data;
 858                         if(btom==&(atom[i]))
 859                                 continue;
 860                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 861                         d=v3_absolute_square(&distance);        /* 1/r^2 */
 862                         if(d<=moldyn->cutoff_square) {
 863                                 d=1.0/d;        /* 1/r^2 */
 864                                 h1=d*d;         /* 1/r^4 */
 865                                 h2*=d;          /* 1/r^6 */
 866                                 h1=h2*h2;       /* 1/r^12 */
 867                                 u+=eps*(sig12*h1-sig6*h2);
 868                                 h2*=d;          /* 1/r^8 */
 869                                 h1*=d;          /* 1/r^14 */
 870                                 h2*=6*sig6;
 871                                 h1*=12*sig12;
 872                                 d=-h1+h2;
 873                                 d*=eps;
 874                                 v3_scale(&force,&distance,d);
 875                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 876                         }
 877                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 878
 879                 /* neighbours not doing boundary condition overflow */
 880                 for(j=1;j<c;j++) {
 881                         this=&(neighbour[j]);
 882                         list_reset(this); /* there might not be a single atom */
 883                         if(this->start!=NULL) {
 884
 885                         do {
 886                                 btom=this->current->data;
 887                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 888                                 d=v3_absolute_square(&distance); /* r^2 */
 889                                 if(d<=moldyn->cutoff_square) {
 890                                         d=1.0/d;        /* 1/r^2 */
 891                                         h1=d*d;         /* 1/r^4 */
 892                                         h2*=d;          /* 1/r^6 */
 893                                         h1=h2*h2;       /* 1/r^12 */
 894                                         u+=eps*(sig12*h1-sig6*h2);
 895                                         h2*=d;          /* 1/r^8 */
 896                                         h1*=d;          /* 1/r^14 */
 897                                         h2*=6*sig6;
 898                                         h1*=12*sig12;
 899                                         d=-h1+h2;
 900                                         d*=eps;
 901                                         v3_scale(&force,&distance,d);
 902                                         v3_add(&(atom[i].f),&(atom[i].f),
 903                                                &force);
 904                                 }
 905                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 906
 907                         }
 908                 }
 909
 910                 /* neighbours due to boundary conditions */
 911                 for(j=c;j<27;j++) {
 912                         this=&(neighbour[j]);
 913                         list_reset(this); /* check boundary conditions */
 914                         if(this->start!=NULL) {
 915
 916                         do {
 917                                 btom=this->current->data;
 918                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 919                                 v3_per_bound(&distance,&(moldyn->dim));
 920                                 d=v3_absolute_square(&distance); /* r^2 */
 921                                 if(d<=moldyn->cutoff_square) {
 922                                         d=1.0/d;        /* 1/r^2 */
 923                                         h1=d*d;         /* 1/r^4 */
 924                                         h2*=d;          /* 1/r^6 */
 925                                         h1=h2*h2;       /* 1/r^12 */
 926                                         u+=eps*(sig12*h1-sig6*h2);
 927                                         h2*=d;          /* 1/r^8 */
 928                                         h1*=d;          /* 1/r^14 */
 929                                         h2*=6*sig6;
 930                                         h1*=12*sig12;
 931                                         d=-h1+h2;
 932                                         d*=eps;
 933                                         v3_scale(&force,&distance,d);
 934                                         v3_add(&(atom[i].f),&(atom[i].f),
 935                                                &force);
 936                                 }
 937                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 938
 939                         }
 940                 }
 941         }
 942
 943         moldyn->energy=0.5*u;
 944
 945         return 0;
 946 }
 947