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