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<nx*ny*nz;i++)
 482                 list_destroy(&(moldyn->lc.subcell[i]));
 483
 484         for(count=0;count<moldyn->count;count++) {
 485                 i=atom[count].r.x/lc->x;
 486                 j=atom[count].r.y/lc->y;
 487                 k=atom[count].r.z/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         cell[0]=lc->subcell[i+j*nx+k*a];
 514         for(ci=-1;ci<=1;ci++) {
 515                 bx=0;
 516                 x=i+ci;
 517                 if((x<0)||(x>=nx)) {
 518                         x=(x+nx)%nx;
 519                         bx=1;
 520                 }
 521                 for(cj=-1;cj<=1;cj++) {
 522                         by=0;
 523                         y=j+cj;
 524                         if((y<0)||(y>=ny)) {
 525                                 y=(y+ny)%ny;
 526                                 by=1;
 527                         }
 528                         for(ck=-1;ck<=1;ck++) {
 529                                 bz=0;
 530                                 z=k+ck;
 531                                 if((z<0)||(z>=nz)) {
 532                                         z=(z+nz)%nz;
 533                                         bz=1;
 534                                 }
 535                                 if(!(x|y|z)) continue;
 536                                 if(bx|by|bz) {
 537                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
 538                                 }
 539                                 else {
 540                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
 541                                 }
 542                         }
 543                 }
 544         }
 545
 546         return count2;
 547 }
 548
 549 int link_cell_shutdown(t_moldyn *moldyn) {
 550
 551         int i;
 552         t_linkcell *lc;
 553
 554         lc=&(moldyn->lc);
 555
 556         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
 557                 list_shutdown(&(moldyn->lc.subcell[i]));
 558
 559         if(lc->listfd) close(lc->listfd);
 560
 561         return 0;
 562 }
 563
 564 /*
 565  *
 566  * 'integration of newtons equation' - algorithms
 567  *
 568  */
 569
 570 /* start the integration */
 571
 572 int moldyn_integrate(t_moldyn *moldyn) {
 573
 574         int i;
 575         unsigned int e,m,s,d,v;
 576         t_3dvec p;
 577
 578         int fd;
 579         char fb[128];
 580
 581         /* logging & visualization */
 582         e=moldyn->ewrite;
 583         m=moldyn->mwrite;
 584         s=moldyn->swrite;
 585         d=moldyn->dwrite;
 586         v=moldyn->vwrite;
 587
 588         if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
 589                 printf("[moldyn] warning, lv system not initialized\n");
 590                 return -1;
 591         }
 592
 593         /* sqaure of some variables */
 594         moldyn->tau_square=moldyn->tau*moldyn->tau;
 595         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 596
 597         /* calculate initial forces */
 598         moldyn->potential_force_function(moldyn);
 599
 600         for(i=0;i<moldyn->time_steps;i++) {
 601                 /* show runs */
 602                 printf(".");
 603
 604                 /* neighbour list update */
 605                 link_cell_update(moldyn);
 606
 607                 /* integration step */
 608                 moldyn->integrate(moldyn);
 609
 610                 /* check for log & visualization */
 611                 if(e) {
 612                         if(!(i%e))
 613                                 dprintf(moldyn->efd,
 614                                         "%.15f %.45f\n",i*moldyn->tau,
 615                                         get_total_energy(moldyn));
 616                 }
 617                 if(m) {
 618                         if(!(i%m)) {
 619                                 p=get_total_p(moldyn->atom,moldyn->count);
 620                                 dprintf(moldyn->mfd,
 621                                         "%.15f %.45f\n",i*moldyn->tau,
 622                                         v3_norm(&p));
 623                         }
 624                 }
 625                 if(s) {
 626                         if(!(i%s)) {
 627                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
 628                                          moldyn->t,i*moldyn->tau);
 629                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
 630                                 if(fd<0) perror("[moldyn] save fd open");
 631                                 else {
 632                                         write(fd,moldyn,sizeof(t_moldyn));
 633                                         write(fd,moldyn->atom,
 634                                               moldyn->count*sizeof(t_atom));
 635                                 }
 636                         }
 637                 }
 638                 if(d) {
 639                         if(!(i%d))
 640                                 write(moldyn->dfd,moldyn->atom,
 641                                       moldyn->count*sizeof(t_atom));
 642
 643                 }
 644                 if(v) {
 645                         if(!(i%v))
 646                                 visual_atoms(moldyn->visual,i*moldyn->tau,
 647                                              moldyn->atom,moldyn->count);
 648                 }
 649         }
 650
 651         return 0;
 652 }
 653
 654 /* velocity verlet */
 655
 656 int velocity_verlet(t_moldyn *moldyn) {
 657
 658         int i,count;
 659         double tau,tau_square;
 660         t_3dvec delta;
 661         t_atom *atom;
 662
 663         atom=moldyn->atom;
 664         count=moldyn->count;
 665         tau=moldyn->tau;
 666         tau_square=moldyn->tau_square;
 667
 668         for(i=0;i<count;i++) {
 669                 /* new positions */
 670                 v3_scale(&delta,&(atom[i].v),tau);
 671                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 672                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 673                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 674                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
 675
 676                 /* velocities */
 677                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 678                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 679         }
 680
 681         /* forces depending on chosen potential */
 682         moldyn->potential_force_function(moldyn);
 683
 684         for(i=0;i<count;i++) {
 685                 /* again velocities */
 686                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 687                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 688         }
 689
 690         return 0;
 691 }
 692
 693
 694 /*
 695  *
 696  * potentials & corresponding forces
 697  *
 698  */
 699
 700 /* harmonic oscillator potential and force */
 701
 702 int harmonic_oscillator(t_moldyn *moldyn) {
 703
 704         t_ho_params *params;
 705         t_atom *atom,*btom;
 706         t_linkcell *lc;
 707         t_list *this,neighbour[27];
 708         int i,j,c;
 709         int count;
 710         t_3dvec force,distance;
 711         double d,u;
 712         double sc,equi_dist;
 713         int ni,nj,nk;
 714
 715         params=moldyn->pot_params;
 716         atom=moldyn->atom;
 717         lc=&(moldyn->lc);
 718         sc=params->spring_constant;
 719         equi_dist=params->equilibrium_distance;
 720         count=moldyn->count;
 721
 722         u=0.0;
 723         for(i=0;i<count;i++) {
 724                 /* determine cell + neighbours */
 725                 ni=atom[i].r.x/lc->x;
 726                 nj=atom[i].r.y/lc->y;
 727                 nk=atom[i].r.z/lc->z;
 728                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 729
 730                 /* processing cell of atom i */
 731                 this=&(neighbour[0]);
 732                 list_reset(this); /* list has 1 element at minimum */
 733                 do {
 734                         btom=this->current->data;
 735                         if(btom==&(atom[i]))
 736                                 continue;
 737                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 738                         d=v3_norm(&distance);
 739                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 740                         v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d)));
 741                         v3_add(&(atom[i].f),&(atom[i].f),&force);
 742                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 743
 744                 /* neighbours not doing boundary condition overflow */
 745                 for(j=1;j<c;j++) {
 746                         this=&(neighbour[j]);
 747                         list_reset(this); /* there might not be a single atom */
 748                         if(this->start!=NULL) {
 749
 750                         do {
 751                                 btom=this->current->data;
 752                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 753                                 d=v3_norm(&distance);
 754                                 if(d<=moldyn->cutoff) {
 755                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 756                                         v3_scale(&force,&distance,
 757                                                  -sc*(1.0-(equi_dist/d)));
 758                                         v3_add(&(atom[i].f),&(atom[i].f),
 759                                                &force);
 760                                 }
 761                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 762
 763                         }
 764                 }
 765
 766                 /* neighbours due to boundary conditions */
 767                 for(j=c;j<27;j++) {
 768                         this=&(neighbour[j]);
 769                         list_reset(this); /* check boundary conditions */
 770                         if(this->start!=NULL) {
 771
 772                         do {
 773                                 btom=this->current->data;
 774                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 775                                 v3_per_bound(&distance,&(moldyn->dim));
 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         moldyn->energy=u;
 791
 792         return 0;
 793 }
 794
 795 /* lennard jones potential & force for one sort of atoms */
 796
 797 int lennard_jones(t_moldyn *moldyn) {
 798
 799         t_lj_params *params;
 800         t_atom *atom,*btom;
 801         t_linkcell *lc;
 802         t_list *this,neighbour[27];
 803         int i,j,c;
 804         int count;
 805         t_3dvec force,distance;
 806         double d,h1,h2,u;
 807         double eps,sig6,sig12;
 808         int ni,nj,nk;
 809
 810         params=moldyn->pot_params;
 811         atom=moldyn->atom;
 812         lc=&(moldyn->lc);
 813         count=moldyn->count;
 814         eps=params->epsilon4;
 815         sig6=params->sigma6;
 816         sig12=params->sigma12;
 817
 818         u=0.0;
 819         for(i=0;i<count;i++) {
 820                 /* determine cell + neighbours */
 821                 ni=atom[i].r.x/lc->x;
 822                 nj=atom[i].r.y/lc->y;
 823                 nk=atom[i].r.z/lc->z;
 824                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
 825
 826                 /* processing cell of atom i */
 827                 this=&(neighbour[0]);
 828                 list_reset(this); /* list has 1 element at minimum */
 829                 do {
 830                         btom=this->current->data;
 831                         if(btom==&(atom[i]))
 832                                 continue;
 833                         v3_sub(&distance,&(atom[i].r),&(btom->r));
 834                         d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
 835                         h1=d*d;         /* 1/r^4 */
 836                         h2*=d;          /* 1/r^6 */
 837                         h1=h2*h2;       /* 1/r^12 */
 838                         u+=eps*(sig12*h1-sig6*h2);
 839                         h2*=d;          /* 1/r^8 */
 840                         h1*=d;          /* 1/r^14 */
 841                         h2*=6*sig6;
 842                         h1*=12*sig12;
 843                         d=-h1+h2;
 844                         d*=eps;
 845                         v3_scale(&force,&distance,d);
 846                         v3_add(&(atom[i].f),&(atom[i].f),&force);
 847                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 848
 849                 /* neighbours not doing boundary condition overflow */
 850                 for(j=1;j<c;j++) {
 851                         this=&(neighbour[j]);
 852                         list_reset(this); /* there might not be a single atom */
 853                         if(this->start!=NULL) {
 854
 855                         do {
 856                                 btom=this->current->data;
 857                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 858                                 d=v3_absolute_square(&distance); /* r^2 */
 859                                 if(d<=moldyn->cutoff_square) {
 860                                         d=1.0/d;        /* 1/r^2 */
 861                                         h1=d*d;         /* 1/r^4 */
 862                                         h2*=d;          /* 1/r^6 */
 863                                         h1=h2*h2;       /* 1/r^12 */
 864                                         u+=eps*(sig12*h1-sig6*h2);
 865                                         h2*=d;          /* 1/r^8 */
 866                                         h1*=d;          /* 1/r^14 */
 867                                         h2*=6*sig6;
 868                                         h1*=12*sig12;
 869                                         d=-h1+h2;
 870                                         d*=eps;
 871                                         v3_scale(&force,&distance,d);
 872                                         v3_add(&(atom[i].f),&(atom[i].f),
 873                                                &force);
 874                                 }
 875                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 876
 877                         }
 878                 }
 879
 880                 /* neighbours due to boundary conditions */
 881                 for(j=c;j<27;j++) {
 882                         this=&(neighbour[j]);
 883                         list_reset(this); /* check boundary conditions */
 884                         if(this->start!=NULL) {
 885
 886                         do {
 887                                 btom=this->current->data;
 888                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
 889                                 v3_per_bound(&distance,&(moldyn->dim));
 890                                 d=v3_absolute_square(&distance); /* r^2 */
 891                                 if(d<=moldyn->cutoff_square) {
 892                                         d=1.0/d;        /* 1/r^2 */
 893                                         h1=d*d;         /* 1/r^4 */
 894                                         h2*=d;          /* 1/r^6 */
 895                                         h1=h2*h2;       /* 1/r^12 */
 896                                         u+=eps*(sig12*h1-sig6*h2);
 897                                         h2*=d;          /* 1/r^8 */
 898                                         h1*=d;          /* 1/r^14 */
 899                                         h2*=6*sig6;
 900                                         h1*=12*sig12;
 901                                         d=-h1+h2;
 902                                         d*=eps;
 903                                         v3_scale(&force,&distance,d);
 904                                         v3_add(&(atom[i].f),&(atom[i].f),
 905                                                &force);
 906                                 }
 907                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 908
 909                         }
 910                 }
 911         }
 912
 913         moldyn->energy=u;
 914
 915         return 0;
 916 }
 917