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("-R <runs> (%d)\n",MOLDYN_RUNS);
  39         printf(" -- integration algo --\n");
  40         printf("  -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
  41         printf("     0: velocity verlet\n");
  42         printf(" -- potential --\n");
  43         printf("  -P <number> <param1 param2 ...>\n");
  44         printf("     0: harmonic oscillator\n");
  45         printf("        param1: spring constant\n");
  46         printf("        param2: equilibrium distance\n");
  47         printf("     1: lennard jones\n");
  48         printf("        param1: epsilon\n");
  49         printf("        param2: sigma\n");
  50         printf("\n");
  51
  52         return 0;
  53 }
  54
  55 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
  56
  57         int i;
  58         t_ho_params hop;
  59         t_lj_params ljp;
  60         double s,e;
  61
  62         memset(moldyn,0,sizeof(t_moldyn));
  63
  64         /* default values */
  65         moldyn->t=MOLDYN_TEMP;
  66         moldyn->tau=MOLDYN_TAU;
  67         moldyn->time_steps=MOLDYN_RUNS;
  68         moldyn->integrate=velocity_verlet;
  69         moldyn->potential_force_function=lennard_jones;
  70
  71         /* parse argv */
  72         for(i=1;i<argc;i++) {
  73                 if(argv[i][0]=='-') {
  74                         switch(argv[i][1]){
  75                                 case 'E':
  76                                         moldyn->ewrite=atoi(argv[++i]);
  77                                         strncpy(moldyn->efb,argv[++i],64);
  78                                         break;
  79                                 case 'M':
  80                                         moldyn->mwrite=atoi(argv[++i]);
  81                                         strncpy(moldyn->mfb,argv[++i],64);
  82                                         break;
  83                                 case 'D':
  84                                         moldyn->dwrite=atoi(argv[++i]);
  85                                         strncpy(moldyn->dfb,argv[++i],64);
  86                                         break;
  87                                 case 'S':
  88                                         moldyn->swrite=atoi(argv[++i]);
  89                                         strncpy(moldyn->sfb,argv[++i],64);
  90                                         break;
  91                                 case 'V':
  92                                         moldyn->vwrite=atoi(argv[++i]);
  93                                         strncpy(moldyn->vfb,argv[++i],64);
  94                                         break;
  95                                 case 'T':
  96                                         moldyn->t=atof(argv[++i]);
  97                                         break;
  98                                 case 't':
  99                                         moldyn->tau=atof(argv[++i]);
 100                                         break;
 101                                 case 'R':
 102                                         moldyn->time_steps=atoi(argv[++i]);
 103                                         break;
 104                                 case 'I':
 105         /* integration algorithm */
 106         switch(atoi(argv[++i])) {
 107                 case MOLDYN_INTEGRATE_VERLET:
 108                         moldyn->integrate=velocity_verlet;
 109                         break;
 110                 default:
 111                         printf("unknown integration algo %s\n",argv[i]);
 112                         moldyn_usage(argv);
 113                         return -1;
 114         }
 115
 116                                 case 'P':
 117         /* potential + params */
 118         switch(atoi(argv[++i])) {
 119                 case MOLDYN_POTENTIAL_HO:
 120                         hop.spring_constant=atof(argv[++i]);
 121                         hop.equilibrium_distance=atof(argv[++i]);
 122                         moldyn->pot_params=malloc(sizeof(t_ho_params));
 123                         memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
 124                         moldyn->potential_force_function=harmonic_oscillator;
 125                         break;
 126                 case MOLDYN_POTENTIAL_LJ:
 127                         e=atof(argv[++i]);
 128                         s=atof(argv[++i]);
 129                         ljp.epsilon4=4*e;
 130                         ljp.sigma6=s*s*s*s*s*s;
 131                         ljp.sigma12=ljp.sigma6*ljp.sigma6;
 132                         moldyn->pot_params=malloc(sizeof(t_lj_params));
 133                         memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
 134                         moldyn->potential_force_function=lennard_jones;
 135                         break;
 136                 default:
 137                         printf("unknown potential %s\n",argv[i]);
 138                         moldyn_usage(argv);
 139                         return -1;
 140         }
 141
 142                                 default:
 143                                         printf("unknown option %s\n",argv[i]);
 144                                         moldyn_usage(argv);
 145                                         return -1;
 146                         }
 147                 } else {
 148                         moldyn_usage(argv);
 149                         return -1;
 150                 }
 151         }
 152
 153         return 0;
 154 }
 155
 156 int moldyn_log_init(t_moldyn *moldyn,void *v) {
 157
 158         moldyn->lvstat=0;
 159         t_visual *vis;
 160
 161         vis=v;
 162
 163         if(moldyn->ewrite) {
 164                 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
 165                 if(moldyn->efd<0) {
 166                         perror("[moldyn] efd open");
 167                         return moldyn->efd;
 168                 }
 169                 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
 170                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
 171         }
 172
 173         if(moldyn->mwrite) {
 174                 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
 175                 if(moldyn->mfd<0) {
 176                         perror("[moldyn] mfd open");
 177                         return moldyn->mfd;
 178                 }
 179                 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
 180                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
 181         }
 182
 183         if(moldyn->swrite)
 184                 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
 185
 186         if(moldyn->dwrite) {
 187                 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
 188                 if(moldyn->dfd<0) {
 189                         perror("[moldyn] dfd open");
 190                         return moldyn->dfd;
 191                 }
 192                 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
 193                 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
 194         }
 195
 196         if((moldyn->vwrite)&&(vis)) {
 197                 moldyn->visual=vis;
 198                 visual_init(vis,moldyn->vfb);
 199                 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
 200         }
 201
 202         moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
 203
 204         return 0;
 205 }
 206
 207 int moldyn_shutdown(t_moldyn *moldyn) {
 208
 209         if(moldyn->efd) close(moldyn->efd);
 210         if(moldyn->mfd) close(moldyn->efd);
 211         if(moldyn->dfd) close(moldyn->efd);
 212         if(moldyn->visual) visual_tini(moldyn->visual);
 213
 214         return 0;
 215 }
 216
 217 int create_lattice(unsigned char type,int element,double mass,double lc,
 218                    int a,int b,int c,t_atom **atom) {
 219
 220         int count;
 221         int ret;
 222         t_3dvec origin;
 223
 224         count=a*b*c;
 225
 226         if(type==FCC) count*=4;
 227         if(type==DIAMOND) count*=8;
 228
 229         *atom=malloc(count*sizeof(t_atom));
 230         if(*atom==NULL) {
 231                 perror("malloc (atoms)");
 232                 return -1;
 233         }
 234
 235         v3_zero(&origin);
 236
 237         switch(type) {
 238                 case FCC:
 239                         ret=fcc_init(a,b,c,lc,*atom,&origin);
 240                         break;
 241                 case DIAMOND:
 242                         ret=diamond_init(a,b,c,lc,*atom,&origin);
 243                         break;
 244                 default:
 245                         printf("unknown lattice type (%02x)\n",type);
 246                         return -1;
 247         }
 248
 249         /* debug */
 250         if(ret!=count) {
 251                 printf("ok, there is something wrong ...\n");
 252                 printf("calculated -> %d atoms\n",count);
 253                 printf("created -> %d atoms\n",ret);
 254                 return -1;
 255         }
 256
 257         while(count) {
 258                 (*atom)[count-1].element=element;
 259                 (*atom)[count-1].mass=mass;
 260                 count-=1;
 261         }
 262
 263         return ret;
 264 }
 265
 266 int destroy_lattice(t_atom *atom) {
 267
 268         if(atom) free(atom);
 269
 270         return 0;
 271 }
 272
 273 int thermal_init(t_moldyn *moldyn,t_random *random) {
 274
 275         /*
 276          * - gaussian distribution of velocities
 277          * - zero total momentum
 278          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 279          */
 280
 281         int i;
 282         double v,sigma;
 283         t_3dvec p_total,delta;
 284         t_atom *atom;
 285
 286         atom=moldyn->atom;
 287
 288         /* gaussian distribution of velocities */
 289         v3_zero(&p_total);
 290         for(i=0;i<moldyn->count;i++) {
 291                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
 292                 /* x direction */
 293                 v=sigma*rand_get_gauss(random);
 294                 atom[i].v.x=v;
 295                 p_total.x+=atom[i].mass*v;
 296                 /* y direction */
 297                 v=sigma*rand_get_gauss(random);
 298                 atom[i].v.y=v;
 299                 p_total.y+=atom[i].mass*v;
 300                 /* z direction */
 301                 v=sigma*rand_get_gauss(random);
 302                 atom[i].v.z=v;
 303                 p_total.z+=atom[i].mass*v;
 304         }
 305
 306         /* zero total momentum */
 307         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 308         for(i=0;i<moldyn->count;i++) {
 309                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 310                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 311         }
 312
 313         /* velocity scaling */
 314         scale_velocity(moldyn);
 315
 316         return 0;
 317 }
 318
 319 int scale_velocity(t_moldyn *moldyn) {
 320
 321         int i;
 322         double e,c;
 323         t_atom *atom;
 324
 325         atom=moldyn->atom;
 326
 327         /*
 328          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 329          */
 330         e=0.0;
 331         for(i=0;i<moldyn->count;i++)
 332                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 333         c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
 334         for(i=0;i<moldyn->count;i++)
 335                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 336
 337         return 0;
 338 }
 339
 340 double get_e_kin(t_atom *atom,int count) {
 341
 342         int i;
 343         double e;
 344
 345         e=0.0;
 346
 347         for(i=0;i<count;i++) {
 348                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 349         }
 350
 351         return e;
 352 }
 353
 354 double get_e_pot(t_moldyn *moldyn) {
 355
 356         return moldyn->energy;
 357 }
 358
 359 double get_total_energy(t_moldyn *moldyn) {
 360
 361         double e;
 362
 363         e=get_e_kin(moldyn->atom,moldyn->count);
 364         e+=get_e_pot(moldyn);
 365
 366         return e;
 367 }
 368
 369 t_3dvec get_total_p(t_atom *atom, int count) {
 370
 371         t_3dvec p,p_total;
 372         int i;
 373
 374         v3_zero(&p_total);
 375         for(i=0;i<count;i++) {
 376                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 377                 v3_add(&p_total,&p_total,&p);
 378         }
 379
 380         return p_total;
 381 }
 382
 383 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
 384
 385         double tau;
 386
 387         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
 388         tau*=1.0E-9;
 389         if(tau<moldyn->tau)
 390                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
 391                        moldyn->tau,tau);
 392
 393         return tau;
 394 }
 395
 396 /*
 397  * numerical tricks
 398  */
 399
 400 /* linked list / cell method */
 401
 402 int link_cell_init(t_moldyn *moldyn) {
 403
 404         t_linkcell *lc;
 405
 406         lc=&(moldyn->lc);
 407
 408         /* partitioning the md cell */
 409         lc->nx=moldyn->dim.x/moldyn->cutoff;
 410         lc->x=moldyn->dim.x/lc->nx;
 411         lc->ny=moldyn->dim.y/moldyn->cutoff;
 412         lc->y=moldyn->dim.y/lc->ny;
 413         lc->nz=moldyn->dim.z/moldyn->cutoff;
 414         lc->z=moldyn->dim.z/lc->nz;
 415
 416         lc->subcell=malloc(lc->nx*lc->ny*lc->nz*sizeof(t_list));
 417
 418         link_cell_update(moldyn);
 419
 420         return 0;
 421 }
 422
 423 int link_cell_update(t_moldyn *moldyn) {
 424
 425         int count,i,j,k;
 426         int nx,ny,nz;
 427         t_atom *atom;
 428         t_linkcell *lc;
 429
 430         atom=moldyn->atom;
 431         lc=&(moldyn->lc);
 432
 433         nx=lc->nx;
 434         ny=lc->ny;
 435         nz=lc->nz;
 436
 437         for(i=0;i<nx*ny*nz;i++)
 438                 list_destroy(&(moldyn->lc.subcell[i]));
 439
 440         for(count=0;count<moldyn->count;count++) {
 441                 i=atom[count].r.x/lc->x;
 442                 j=atom[count].r.y/lc->y;
 443                 k=atom[count].r.z/lc->z;
 444                 list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
 445                                        &(atom[count]));
 446         }
 447
 448         return 0;
 449 }
 450
 451 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 452
 453         t_linkcell *lc;
 454         int a;
 455         int count1,count2;
 456         int ci,cj,ck;
 457         int nx,ny,nz;
 458         int x,y,z;
 459         unsigned char bx,by,bz;
 460
 461         lc=&(moldyn->lc);
 462         nx=lc->nx;
 463         ny=lc->ny;
 464         nx=lc->nz;
 465         count1=1;
 466         count2=27;
 467         a=nx*ny;
 468
 469         cell[0]=lc->subcell[i+j*nx+k*a];
 470         for(ci=-1;ci<=1;ci++) {
 471                 bx=0;
 472                 x=i+ci;
 473                 if((x<0)||(x>=nx)) {
 474                         x=(x+nx)%nx;
 475                         bx=1;
 476                 }
 477                 for(cj=-1;cj<=1;cj++) {
 478                         by=0;
 479                         y=j+cj;
 480                         if((y<0)||(y>=ny)) {
 481                                 y=(y+ny)%ny;
 482                                 by=1;
 483                         }
 484                         for(ck=-1;ck<=1;ck++) {
 485                                 bz=0;
 486                                 z=k+ck;
 487                                 if((z<0)||(z>=nz)) {
 488                                         z=(z+nz)%nz;
 489                                         bz=1;
 490                                 }
 491                                 if(!(x|y|z)) continue;
 492                                 if(bx|by|bz) {
 493                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
 494                                 }
 495                                 else {
 496                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
 497                                 }
 498                         }
 499                 }
 500         }
 501
 502         return count2;
 503 }
 504
 505 int link_cell_shutdown(t_moldyn *moldyn) {
 506
 507         int i;
 508         t_linkcell *lc;
 509
 510         lc=&(moldyn->lc);
 511
 512         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
 513                 list_shutdown(&(moldyn->lc.subcell[i]));
 514
 515         return 0;
 516 }
 517
 518 /*
 519  *
 520  * 'integration of newtons equation' - algorithms
 521  *
 522  */
 523
 524 /* start the integration */
 525
 526 int moldyn_integrate(t_moldyn *moldyn) {
 527
 528         int i;
 529         unsigned int e,m,s,d,v;
 530         t_3dvec p;
 531
 532         int fd;
 533         char fb[128];
 534
 535         /* logging & visualization */
 536         e=moldyn->ewrite;
 537         m=moldyn->mwrite;
 538         s=moldyn->swrite;
 539         d=moldyn->dwrite;
 540         v=moldyn->vwrite;
 541
 542         if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
 543                 printf("[moldyn] warning, lv system not initialized\n");
 544                 return -1;
 545         }
 546
 547         /* sqaure of some variables */
 548         moldyn->tau_square=moldyn->tau*moldyn->tau;
 549         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 550
 551         /* create the neighbour list */
 552         link_cell_update(moldyn);
 553
 554         /* calculate initial forces */
 555         moldyn->potential_force_function(moldyn);
 556
 557         for(i=0;i<moldyn->time_steps;i++) {
 558                 /* show runs */
 559                 printf(".");
 560
 561                 /* neighbour list update */
 562                 link_cell_update(moldyn);
 563
 564                 /* integration step */
 565                 moldyn->integrate(moldyn);
 566
 567                 /* check for log & visualization */
 568                 if(e) {
 569                         if(!(i%e))
 570                                 dprintf(moldyn->efd,
 571                                         "%.15f %.45f\n",i*moldyn->tau,
 572                                         get_total_energy(moldyn));
 573                 }
 574                 if(m) {
 575                         if(!(i%m)) {
 576                                 p=get_total_p(moldyn->atom,moldyn->count);
 577                                 dprintf(moldyn->mfd,
 578                                         "%.15f %.45f\n",i*moldyn->tau,
 579                                         v3_norm(&p));
 580                         }
 581                 }
 582                 if(s) {
 583                         if(!(i%s)) {
 584                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
 585                                          moldyn->t,i*moldyn->tau);
 586                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
 587                                 if(fd<0) perror("[moldyn] save fd open");
 588                                 else {
 589                                         write(fd,moldyn,sizeof(t_moldyn));
 590                                         write(fd,moldyn->atom,
 591                                               moldyn->count*sizeof(t_atom));
 592                                 }
 593                         }
 594                 }
 595                 if(d) {
 596                         if(!(i%d))
 597                                 write(moldyn->dfd,moldyn->atom,
 598                                       moldyn->count*sizeof(t_atom));
 599
 600                 }
 601                 if(v) {
 602                         if(!(i%v))
 603                                 visual_atoms(moldyn->visual,i*moldyn->tau,
 604                                              moldyn->atom,moldyn->count);
 605                 }
 606         }
 607
 608         return 0;
 609 }
 610
 611 /* velocity verlet */
 612
 613 int velocity_verlet(t_moldyn *moldyn) {
 614
 615         int i,count;
 616         double tau,tau_square;
 617         t_3dvec delta;
 618         t_atom *atom;
 619
 620         atom=moldyn->atom;
 621         count=moldyn->count;
 622         tau=moldyn->tau;
 623         tau_square=moldyn->tau_square;
 624
 625         for(i=0;i<count;i++) {
 626                 /* new positions */
 627                 v3_scale(&delta,&(atom[i].v),tau);
 628                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 629                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 630                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 631                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
 632
 633                 /* velocities */
 634                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 635                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 636         }
 637
 638         /* forces depending on chosen potential */
 639         moldyn->potential_force_function(moldyn);
 640
 641         for(i=0;i<count;i++) {
 642                 /* again velocities */
 643                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 644                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 645         }
 646
 647         return 0;
 648 }
 649
 650
 651 /*
 652  *
 653  * potentials & corresponding forces
 654  *
 655  */
 656
 657 /* harmonic oscillator potential and force */
 658
 659 int harmonic_oscillator(t_moldyn *moldyn) {
 660
 661         t_ho_params *params;
 662         t_atom *atom,*btom;
 663         t_linkcell *lc;
 664         t_list *this,neighbour[27];
 665         int i,j,c;
 666         int count;
 667         t_3dvec force,distance;
 668         double d,u;
 669         double sc,equi_dist;
 670         int ni,nj,nk;
 671
 672         params=moldyn->pot_params;
 673         atom=moldyn->atom;
 674         lc=&(moldyn->lc);
 675         sc=params->spring_constant;
 676         equi_dist=params->equilibrium_distance;
 677         count=moldyn->count;
 678
 679         u=0.0;
 680         for(i=0;i<count;i++) {
 681                 /* determine cell + neighbours */
 682                 ni=atom[i].r.x/lc->x;
 683                 nj=atom[i].r.y/lc->y;
 684                 nk=atom[i].r.z/lc->z;
 685                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 686
 687                 /* processing cell of atom i */
 688                 this=&(neighbour[0]);
 689                 list_reset(this); /* list has 1 element at minimum */
 690                 do {
 691                         btom=this->current->data;
 692                         if(btom==&(atom[i]))
 693                                 continue;
 694                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 695                         d=v3_norm(&distance);
 696                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 697                         v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d)));
 698                         v3_add(&(atom[i].f),&(atom[i].f),&force);
 699                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 700
 701                 /* neighbours not doing boundary condition overflow */
 702                 for(j=1;j<c;j++) {
 703                         this=&(neighbour[j]);
 704                         list_reset(this); /* there might not be a single atom */
 705                         if(this->start!=NULL) {
 706
 707                         do {
 708                                 btom=this->current->data;
 709                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 710                                 d=v3_norm(&distance);
 711                                 if(d<=moldyn->cutoff) {
 712                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 713                                         v3_scale(&force,&distance,
 714                                                  -sc*(1.0-(equi_dist/d)));
 715                                         v3_add(&(atom[i].f),&(atom[i].f),
 716                                                &force);
 717                                 }
 718                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 719
 720                         }
 721                 }
 722
 723                 /* neighbours due to boundary conditions */
 724                 for(j=c;j<27;j++) {
 725                         this=&(neighbour[j]);
 726                         list_reset(this); /* check boundary conditions */
 727                         if(this->start!=NULL) {
 728
 729                         do {
 730                                 btom=this->current->data;
 731                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 732                                 v3_per_bound(&distance,&(moldyn->dim));
 733                                 d=v3_norm(&distance);
 734                                 if(d<=moldyn->cutoff) {
 735                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 736                                         v3_scale(&force,&distance,
 737                                                  -sc*(1.0-(equi_dist/d)));
 738                                         v3_add(&(atom[i].f),&(atom[i].f),
 739                                                &force);
 740                                 }
 741                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 742
 743                         }
 744                 }
 745         }
 746
 747         moldyn->energy=u;
 748
 749         return 0;
 750 }
 751
 752 /* lennard jones potential & force for one sort of atoms */
 753
 754 int lennard_jones(t_moldyn *moldyn) {
 755
 756         t_lj_params *params;
 757         t_atom *atom,*btom;
 758         t_linkcell *lc;
 759         t_list *this,neighbour[27];
 760         int i,j,c;
 761         int count;
 762         t_3dvec force,distance;
 763         double d,h1,h2,u;
 764         double eps,sig6,sig12;
 765         int ni,nj,nk;
 766
 767         params=moldyn->pot_params;
 768         atom=moldyn->atom;
 769         lc=&(moldyn->lc);
 770         count=moldyn->count;
 771         eps=params->epsilon4;
 772         sig6=params->sigma6;
 773         sig12=params->sigma12;
 774
 775         u=0.0;
 776         for(i=0;i<count;i++) {
 777                 /* determine cell + neighbours */
 778                 ni=atom[i].r.x/lc->x;
 779                 nj=atom[i].r.y/lc->y;
 780                 nk=atom[i].r.z/lc->z;
 781                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 782
 783                 /* processing cell of atom i */
 784                 this=&(neighbour[0]);
 785                 list_reset(this); /* list has 1 element at minimum */
 786                 do {
 787                         btom=this->current->data;
 788                         if(btom==&(atom[i]))
 789                                 continue;
 790                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 791                         d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
 792                         h1=d*d;         /* 1/r^4 */
 793                         h2*=d;          /* 1/r^6 */
 794                         h1=h2*h2;       /* 1/r^12 */
 795                         u+=eps*(sig12*h1-sig6*h2);
 796                         h2*=d;          /* 1/r^8 */
 797                         h1*=d;          /* 1/r^14 */
 798                         h2*=6*sig6;
 799                         h1*=12*sig12;
 800                         d=-h1+h2;
 801                         d*=eps;
 802                         v3_scale(&force,&distance,d);
 803                         v3_add(&(atom[i].f),&(atom[i].f),&force);
 804                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 805
 806                 /* neighbours not doing boundary condition overflow */
 807                 for(j=1;j<c;j++) {
 808                         this=&(neighbour[j]);
 809                         list_reset(this); /* there might not be a single atom */
 810                         if(this->start!=NULL) {
 811
 812                         do {
 813                                 btom=this->current->data;
 814                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 815                                 d=v3_absolute_square(&distance); /* r^2 */
 816                                 if(d<=moldyn->cutoff_square) {
 817                                         d=1.0/d;        /* 1/r^2 */
 818                                         h1=d*d;         /* 1/r^4 */
 819                                         h2*=d;          /* 1/r^6 */
 820                                         h1=h2*h2;       /* 1/r^12 */
 821                                         u+=eps*(sig12*h1-sig6*h2);
 822                                         h2*=d;          /* 1/r^8 */
 823                                         h1*=d;          /* 1/r^14 */
 824                                         h2*=6*sig6;
 825                                         h1*=12*sig12;
 826                                         d=-h1+h2;
 827                                         d*=eps;
 828                                         v3_scale(&force,&distance,d);
 829                                         v3_add(&(atom[i].f),&(atom[i].f),
 830                                                &force);
 831                                 }
 832                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 833
 834                         }
 835                 }
 836
 837                 /* neighbours due to boundary conditions */
 838                 for(j=c;j<27;j++) {
 839                         this=&(neighbour[j]);
 840                         list_reset(this); /* check boundary conditions */
 841                         if(this->start!=NULL) {
 842
 843                         do {
 844                                 btom=this->current->data;
 845                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 846                                 v3_per_bound(&distance,&(moldyn->dim));
 847                                 d=v3_absolute_square(&distance); /* r^2 */
 848                                 if(d<=moldyn->cutoff_square) {
 849                                         d=1.0/d;        /* 1/r^2 */
 850                                         h1=d*d;         /* 1/r^4 */
 851                                         h2*=d;          /* 1/r^6 */
 852                                         h1=h2*h2;       /* 1/r^12 */
 853                                         u+=eps*(sig12*h1-sig6*h2);
 854                                         h2*=d;          /* 1/r^8 */
 855                                         h1*=d;          /* 1/r^14 */
 856                                         h2*=6*sig6;
 857                                         h1*=12*sig12;
 858                                         d=-h1+h2;
 859                                         d*=eps;
 860                                         v3_scale(&force,&distance,d);
 861                                         v3_add(&(atom[i].f),&(atom[i].f),
 862                                                &force);
 863                                 }
 864                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 865
 866                         }
 867                 }
 868         }
 869
 870         moldyn->energy=u;
 871
 872         return 0;
 873 }
 874