//
// Simple MD program made for demonstration purposes
// 
//
// Internal units:
//
// All lengths in Å = 1e-10 m
// All times in fs = 1e-15 s
// All velocities and accelerations in Å/fs and Å/fs^2
// All forces in eV/Å
// All energies in eV
//

// Note that since this code is converted from Fortran, array
// indices run from 1 to N, not from 0 to N-1 as in usually in C

#define MAIN 
#include "global.h"

main()
{  

  // Atom data containing arrays
  // 
  // x, v and a are the atom positions, velocities and accelerations
  // Epot and Ekin are the potential and kinetic energy
  // atomname is simply the atom name
  //
  struct vector   x[MAXAT+1],v[MAXAT+1],a[MAXAT+1];
  double   Epot[MAXAT+1],Ekin[MAXAT+1];
  char atomname[MAXAT+1][5];
  int      neighbourlist[MAXAT*MAXNEIGHBOURS+1];
  
  // Other parameters and variables 
  //
  // time is the current time, tmax the maximum time, deltat the time step
  //   and itime the integer number of the current time step
  // T is the Temperature (in K) and m is the atom mass
  // Ekinsum, Epotsum and Etot are the total energies for all atoms
  // N is the number of atoms
  // nneighbourlistupdate determines the interval between neighbour list
  //    updates, and nmovieoutput the interval between atom position output.
  // rpotcut and rskincut are the potential and neighbour list cutoff radii
  // size is the simulation cell size
  // periodic determines whether periodic boundaries will be used (1=yes, 0=no)
  //
  double   time,tmax,deltat,T,Ekinsum,Epotsum,Etot,m;
  int        N,nneighbourlistupdate,itime,nmovieoutput;
  double   rpotcut,rskincut;
  struct vector   size;
  struct ivector  periodic;

  // Morse potential parameters
  double  morseDe,morsealpha,morseRe;

  // Parameters for Berendsen temperature control [J. Chem. Phys.
  // 81 (1984) 3684]
  // btctau is the Berendsen temperature control tau, initialT the initial
  // temperature and desiredT the desired temperature.
  double btctau,initialT,desiredT;
  
  // Parameters for Berendsen pressure control [J. Chem. Phys.
  //  81 (1984) 3684]
  // bpctau is the Berendsen pressure control tau, bpcbeta 1/B,
  // P the instantaneous pressure, and desiredP the desired pressure
  double bpctau,bpcbeta,P,desiredP,virial,mu;


  int  i;
  
  printf("\n--------------- mdmorse03 V1.2 --------------------\n\n");
  //
  // Initialize the simulation
  //
  ReadParams(&initialT,&m,&size,&tmax,&periodic,  
	     &deltat,&nneighbourlistupdate,&rpotcut,&rskincut,  
	     &morseDe,&morsealpha,&morseRe,&nmovieoutput,
	     &btctau,&T,&bpctau,&bpcbeta,&desiredP);
   
  // btc: Store desired T for temperature control
  desiredT=T;

  ReadAtoms(x,atomname,&N);
  
  // Factor 2 because of energy equipartition theorem
  SetTemperature(v,N,m,2.0*initialT);
  // Check that we got it right
  GetTemperature(v,N,m,&T);
  printf("Initial atom temperature is %g\n",T);

  // Initialize accelerations to 0
  for (i=1;i<=N;i++) {
     a[i].x=0.0; a[i].y=0.0; a[i].z=0.0;
  }
  
  // Start main loop over times
  time=0.0;
  itime=0;
  while (1) {
    if (itime%nneighbourlistupdate==0) {
      UpdateNeighbourlist(x,N,size,periodic,neighbourlist,rskincut);
    }
     
    Solve1(x,v,a,N,size,periodic,deltat);

    GetForces(x,a,N,size,periodic,m,neighbourlist,rpotcut, 
	      morseDe,morsealpha,morseRe,Epot,&virial);

    GetTemperature(v,N,m,&T);

    Solve2(v,a,N,deltat,btctau,T,desiredT);

    // Update time
    time=time+deltat;
     
    // Actual solution done, now compute result quantities

    GetEnergies(v,N,m,&Ekinsum,&Epotsum,&Etot,Ekin,Epot);

    // Compute pressure in units of kbar
    P=(N*kB*T+1.0/3.0*virial*e)/(size.x*size.y*size.z*1e-30)/1e8;

    printf("%s %10.3f %10.3f %12.5f %12.5f %12.5f %12.5f\n","ec ",
	   time,T,Ekinsum/N,Epotsum/N,Etot/N,P);

    if (itime%nmovieoutput==0) {
      printf("%s %10.3f\n","Outputting atom movie at t = ",time);
      WriteAtoms(x,atomname,Ekin,Epot,N,time,size);
    }

    // Do Berendsen pressure control
    if (bpctau > 0.0) {
      mu=pow((1.0-bpcbeta*deltat/bpctau*(desiredP-P)),(1.0/3.0));

      size.x=size.x*mu;
      size.y=size.y*mu;
      size.z=size.z*mu;
      for (i=1;i<=N;i++) {
	x[i].x=x[i].x*mu;
	x[i].y=x[i].y*mu;
	x[i].z=x[i].z*mu;
      }
      if (itime%10==1) {
	printf("bpc %10.6f %10.6f %10.6f %12.5f %12.5f %12.5f\n",
	       size.x,size.y,size.z,size.x*size.y*size.z,P,mu);
      }
    }

    // Check whether we are done
    if (time >= tmax) break;

    itime=itime+1;

  } // End of loop over time steps

  // Always output final energy and atom coordinates
  // using high precision 
  printf("%s %12.5f %12.5f %14.7f %14.7f %14.7f\n","ec ",
	 time,T,Ekinsum/N,Epotsum/N,Etot/N);

  WriteAtoms(x,atomname,Ekin,Epot,N,time,size);

}