% /**************************************************************************************************************************************************
% LJ_2.cpp
% 
% For the two atom unit cell
% 
% Examines the behaviour of an arbitrary (?) number of particles of a Lennard-Jones material
% with periodic boundary conditions. Simulations can be run in either an NVE, NVT or NPT ensemble. Temperature
% is controlled with a Nose-Hoover thermostat. Pressure is controlled in a secret way. The heat current can be calculated in order to predict
% thermal conductivity.
% 
% Can do different types of non-equilibrium simulations to look at model interactions
% 
% Written by Alan McGaughey
% 
% Developed fall 2000
% Last modified (significantly) March 2, 2004
% ***************************************************************************************************************************************************
% Running the program.
fprintf('Running the program. Please be patient.\n');

% Opening all input and output file pointers.
inFileInitial = fopen('initialdata_4.txt', 'r');      % Input files. Read only.
inFilePhi = fopen('phi4.txt', 'r');
inFileA = fopen('a4.txt', 'r');
posFile = fopen('position.xyz', 'w');                 % Output files. Read and write capable.
velFile = fopen('velocity.xyz', 'w');
totEnergyFile = fopen('energy.txt', 'w');
heatCurrentFile = fopen('heatcurrent.txt', 'w');

% Declaring and initializing some global variables.
rstar = 3.4E-10;     % (m)  these are the scales for argon
estar = 1.67E-21;    % (J)
mstar = 6.63E-26;    % (kg)
kB = 1.3806E-23;
N = 32;              % No of particles.
L = zeros(1,3);
global energy_cutoff force_cutoff a

initial_data = 1;    % Treat as bool function. 1 = true. 0 = false.
t= 0;
tstep = 0;
tstep_total = 10000;     % Original = 120000000
tramp = 10;
tstart = 0;
lennard = 5;
delt = 0.002;

m = zeros(1,N);                % Particle masses.
partiden = zeros(1,N);		   % Particle identification for CHIME
deltaR = zeros(1,3);	
global r dF dF2

% Thermostat parameters
Tset = (kB/estar)*0.;	       % Set temperature
tauT = 0.05;	               % Thermostat time constant
s = 1.0;          		       % Dynamical varible used to check thermostat
eta = 0;	                   % Damping factor
Q = 3*Tset*tauT*tauT/2;		   % Used to check thermostat
global s1         			   % Storage of s for energy calculations
	
% Barostat parameters
pressure_on = 0;                         % Treat as bool function. 1 = true. 0 = false.
Pset = (rstar*rstar*rstar)/estar*0;      % Set pressure - number after * is in Pascals
tauP = 0.5;
epsilon = zeros(1,3);
epsilon_bar = zeros(1,3);
pressure = 0;
global volume

% Set the initial values of the coordinates and their momenta.
% xn, pn, pn1, F, xo = position, momentum and force arrays.
xn = zeros(N,3); 
pn = zeros(N,3); 
pn1 = zeros(N,3);
F = zeros(N,3);
xo = zeros(N,3);
Frms = 0;
Frmstot = 0;
Lo2 = zeros(1,3);
global cube_size l_c l_c_2   % Size of cube in terms of particle spacing.
                         
nx = 2;
ny = nx;
nz = nx;           % N = 4*nx*ny*nz must be satisfied

% Place the particles in a cubic lattice
lattice = 'f';     % f = face cenetered cubic
m1 = 1;
m2 = 0.1;

% Setting up the initial locations of the particle position and their type.
if(lattice == 'f')
  l_c = 1.56;
  for i = 0:nx-1
    for j = 0:ny-1
      for k = 0:nz-1
        xn(4*ny*nz*i+4*nz*j+4*k+1, 1) = i*l_c;
        xn(4*ny*nz*i+4*nz*j+4*k+1, 2) = j*l_c;
        xn(4*ny*nz*i+4*nz*j+4*k+1, 3) = k*l_c;
        partiden(1, 4*ny*nz*i+4*nz*j+4*k+1) = 'C';
        m(1, 4*ny*nz*i+4*nz*j+4*k+1) = m1;
        
        xn(4*ny*nz*i+4*nz*j+4*k+2, 1) = i*l_c + l_c/2;
        xn(4*ny*nz*i+4*nz*j+4*k+2, 2) = j*l_c + l_c/2;
		xn(4*ny*nz*i+4*nz*j+4*k+2, 3) = k*l_c;
		partiden(1, 4*ny*nz*i+4*nz*j+4*k+2) = 'S';
		m(1, 4*ny*nz*i+4*nz*j+4*k+2) = m2;

        xn(4*ny*nz*i+4*nz*j+4*k+3, 1) = i*l_c + l_c/2;
		xn(4*ny*nz*i+4*nz*j+4*k+3, 2) = j*l_c;
		xn(4*ny*nz*i+4*nz*j+4*k+3, 3) = k*l_c + l_c/2;
		partiden(1, 4*ny*nz*i+4*nz*j+4*k+3) = 'C';
		m(1, 4*ny*nz*i+4*nz*j+4*k+3) = m1;

        xn(4*ny*nz*i+4*nz*j+4*k+4, 1) = i*l_c;
		xn(4*ny*nz*i+4*nz*j+4*k+4, 2) = j*l_c + l_c/2;
		xn(4*ny*nz*i+4*nz*j+4*k+4, 3) = k*l_c + l_c/2;
		partiden(1, 4*ny*nz*i+4*nz*j+4*k+4) = 'S';
		m(1, 4*ny*nz*i+4*nz*j+4*k+4) = m2;
      end;
    end;
  end; 
end;


% Setting up the cutoff values.
a = 2.5;      % Cutoff radius
energy_cutoff = -4/(a^12) + 4/(a^6);
force_cutoff = 24/(a^7) - 48/(a^13);

% Set the particle momenta (random)
scale = 2E-10;          % This is really small - good for making sure that equilibrium positions are exactly the intial positions
						% starting a run in NVT will get the velocities up to the desired value quickly
mom = zeros(1,3);
rn = 0;

% The reason for the weird rand function is to match it with the C++ rand() function.
for i = 1:N
  for j = 1:3
    rn = ceil(rand*32767);
    if (rn <= 32767/2) pn(i,j) = scale*rn;
    else pn(i,j) = -1*scale*rn;
    end;
  end;
end;

% Read in initialdata_4.txt.
if (initial_data == 1)
  for i = 1:N
    testA = fscanf(inFileInitial, '%c %e %e %e %e %e %e %e');
  end;
end;

% Force the total momentum to be zero in each direction.
for i = 1:N
  for j = 1:3 mom(1,j) = mom(1,j) + pn(i,j); end;
end;

for i = 1:3 mom(1,i) = mom(1,i)/N; end;

for i = 1:N
  for j = 1:3 pn(i,j) = pn(i,j) - mom(1,j); end;
end;

for i = 1:3 mom(1,i) = 0; end;

for i = 1:N
  for j = 1:3 mom(1,j) = mom(1,j) + pn(i,j); end;
end;

% Write the initial particle positions to an output files.
% Write down also the number of data sets for the energy and heatcurrent file.
position_xyz(t,N, partiden, xn, posFile);
fprintf(totEnergyFile, '%d\n', tstep_total/tramp);
fprintf(heatCurrentFile, '%d\n', tstep_total/tramp);

quench = 0;             % Treat as bool function. 1 = true. 0 = false.
Tset1 = (kB/estar)*300.;
Tset2 = (kB/estar)*50.;
time1 = 10000;
ramp = 80000;

TempSet = 1;            % Boolean used to switch from NVT to NVE
block = 1500000;        % Length of block in time steps
block_NVT = 300000;     % Length of intial NVT block
block_count = -1;       % Current block number
global shortrange shortrange_set Tint Tact  % Actual short range energy and desired value. Tact = actual temperature.

aphi = zeros(8,3,2);

for i = 1:8
  aphi(i,1,1) = fscanf(inFileA, '%e', 1);
  aphi(i,1,2) = fscanf(inFilePhi, '%e', 1);
end;

Tset = 0.*(kB/estar);
Tint = 10;

% Scale positions to temp
for i = 1:N
  for j = 1:3
    xo(i,j) = xn(i,j)*aphi(1,1,1)/1.56;
    xn(i,j) = xo(i,j);  % Initial positions are the same as the equilibrium positions.
  end;
end;

L(1,1) = nx*aphi(1, 1, 1);
L(1,2) = L(1,1);
L(1,3) = L(1,1);
Tint = 0;

% Intergate the equations of motion using a Verlet leapfrog algorithm.
while (tstep < tstep_total)
  if(mod(tstep,(10*block)) == 0) 
    Tset = Tset + 10.*kB/estar;
    Tint = Tint + 10;
  end;

  if (mod(tstep,block) == 0)
    block_count = block_count + 1;
    cube_size = nx;
    L(1,1) = nx*aphi(((Tint-10)/10)+1,(nx-2)+1, 0+1);
    L(1,2) = L(1,1);
    L(1,3) = L(1,1);
    shortrange_set = N*aphi(((Tint-10)/10)+1, (nx-2)+1, 1+1);

    Lo2(1,1) = L(1,1)/2;
    Lo2(1,2) = L(1,2)/2;
    Lo2(1,3) = L(1,3)/2;
    l_c = L(1,1)/cube_size;
    l_c_2 = l_c/2.;
    a = 2.5;
    energy_cutoff = -4/a^12 + 4/a^6;
    force_cutoff = 24/a^7 - 48/a^13;
    volume = L(1,1)*L(1,2)*L(1,3);	
  end;	

  % Randomizing velocities.
  if (mod(tstep,block) == 0)	
    for i = 1:N
      for j = 1:3
		rn = ceil(rand*32767);
		if(rn <= 32767/2) pn(i,j) = scale*rn;
		else pn(i,j) = -1.*scale*rn;
		end;
	  end;
	end;		

    for i = 1:3 mom(1,i) = 0; end;

    for i = 1:N
      for j = 1:3 mom(1,j) = mom(1,j) + pn(i,j); end;
    end;
		
    for i = 1:3 mom(1,i) = mom(1,i)/N; end;
   
    for i = 1:N 
	  for j = 1:3 pn(i,j) = pn(i,j) - mom(1,j); end;
    end;
  end;

  % Set the temperature close to the desired value
  KE = KEnergy(pn,m,N);
  if (tstep < block*block_count + block_NVT && tstep >= block*block_count) eta = eta + delt/(tauT*tauT)*(KE/(1.5*(N-1)*Tset)-1.); end; % NVT
  % Start temperature setting procedure, set to true to deactivate (also uncomment else line above)
  if (tstep == block*block_count + block_NVT) TempSet = 0; end;

  if (tstep < block*(block_count+1) && tstep > block*block_count + block_NVT && TempSet == 0)
    shortrange = 0;
    for i = 1:N
      for j = i+1:N
        for k = 1:3
          deltaR(1,k) = xn(i,k) - xn(j,k);
          if (deltaR(1,k) > Lo2(1,k)) deltaR(1,k) = deltaR(1,k) - L(1,k); end;      % Periodic boundary conditions.
          if (deltaR(1,k) < (-1*Lo2(1,k))) deltaR(1,k) = deltaR(1,k) + L(1,k); end;
        end;
  
        if ((deltaR(1,1) < a) && (deltaR(1,2) < a) && (deltaR(1,3) < a)) r = partsep(deltaR);
        else r = 10*a; end;
      
        LJE = LJEnergy(r, energy_cutoff);
        if (r < a) shortrange = shortrange + LJE; end;
      end;
    end;
   
    if (abs(shortrange - shortrange_set) <= 0.001)
      KE = KEnergy(pn,m,N);
      Tact = 2*KE/(3*(N-1));

      for i = 1:N
        for j = 1:3 pn(i,j) = pn(i,j)*sqrt(Tset/Tact); end;     % Scale momenta.
      end;
      
      TempSet = 1;      % Treat as bool function. 1 = true. 0 = false.
      eta = 0;          % Switch to NVE.
    else 
      KE = KEnergy(pn,m,N);
      eta = eta + delt/(tauT*tauT)*(KE/(1.5*(N-1)*Tset)-1); % Stay in NVT if tolerance not met.
    end;
  end;

  t = t + delt;
  tstep = tstep + 1;

  if (quench == 1)
    if (tstep < time1) Tset = Tset1; end;
    if (tstep < time1 + ramp && tstep >= time1) Tset = Tset + (Tset2-Tset1)/ramp; end;
    if (tstep >= time1+ramp) Tset = Tset2; end;
  end;
	
  if (pressure_on == 0) 
	for i = 1:3 
	  epsilon(1,i) = 0;
	  epsilon_bar(1,i) = 0;
    end;
  end;		
			
  % Find the new positions
  for i = 1:N
    for j = 1:3
      xn(i,j) = xn(i,j)*(1+epsilon(1,j)*delt/2)/(1-epsilon(1,j)*delt/2) + (pn(i,j)*delt/m(1,i))/(1-epsilon(1,j)*delt/2);
      if (xn(i,j) < (-1*l_c_2)) xn(i,j) = xn(i,j) + L(1,j); end;     % Periodic boundary conditions
      if (xn(i,j) > (L(1,j)-l_c_2)) xn(i,j) = xn(i,j) - L(1,j); end;
    end; 
  end;		

  % Recalculate contants if barostat in use.
  if (pressure_on == 1) 
    for i = 1:3 
      L(1,i) = L(1,i)*(1+epsilon(1,i)*delt/2)/(1-epsilon(1,i)*delt/2);
      Lo2(1,i) = L(1,i)/2;
	end;

    l_c = L(1,1)/cube_size;
    l_c_2 = l_c/2.;
    a = 2.5;
    energy_cutoff = -4/a^12 + 4/a^6;
    force_cutoff = 24/a^7 - 48/a^13;
    volume = L(1,1)*L(1,2)*L(1,3);	
  end;

  % Set the forces and pressure to zero
  F = zeros(N,3);
  pressure = 0;
  
  % Calculate the forces.
  for i = 1:N
    for j = i+1:N
      for k = 1:3
        deltaR(1,k) = xn(i,k) - xn(j,k);
        if (deltaR(1,k) > Lo2(1,k)) deltaR(1,k) = deltaR(1,k) - L(1,k); end;     % Periodic boundary condtins
        if (deltaR(1,k) < (-1*Lo2(1,k))) deltaR(1,k) = deltaR(1,k) + L(1,k); end; 
      end;
     
      if ((deltaR(1,1) < a) && (deltaR(1,2) < a) && (deltaR(1,3) < a)) r = partsep(deltaR);
      else r = 10*a;    % Otherwise, give 'r' a value greater than the cutoff
      end;   
      
      if(r < a)
        LJF = LJForce(r);
        dF = LJF;
        for k = 1:3
          dF2 = deltaR(1,k)*dF;
          F(i,k) = F(i,k) + dF2;
          F(j,k) = F(j,k) - dF2;
        end;
        pressure = pressure + r*r*dF;  % dF already contains a 1/r term
      end;
    end;
  end;

  KE = KEnergy(pn,m,N);
  pressure = pressure/(3*volume)+ N*2*KE/(3*(N-1))/volume;
  
  % Calculate epsilon if necessary.
  if (pressure_on == 1)
    for i = 1:3
      epsilon_bar(1,i) = epsilon(1,i);
      % Not quite right as it should be the component of the pressure, not the total value
      epsilon(1,i) = epsilon(1,i) + delt/(tauP*tauP)*volume/(N*Tset)*(pressure-Pset); 
      epsilon_bar(1,i) = (epsilon(1,i)+epsilon_bar(1,i))/2;
    end;
  end;

  % Step the momenta.
  % For output calculations.
  if (mod(tstep,lennard) == 0) 
    s1 = s*exp(eta*delt);
    s = (s+s1)/2;

    for i = 1:N 
      for j = 1:3
        pn1(i,j) = pn(i,j)*(1-(eta+epsilon_bar(1,j))*delt/2)/(1+(eta+epsilon_bar(1,j))*delt/2) + F(i,j)*delt/(1+(eta+epsilon_bar(1,j))*delt/2);
        pn(i,j) = (pn(i,j)+pn1(i,j))/2;
	  end;
    end;

    % Update momenta for next loop			
    for i = 1:N
      for j = 1:3 pn(i,j) = pn1(i,j); end;
	end;
    s = s1;
  % Otherwise (normal loop - don't need to average momenta and s - no calculations are being made.
  else
    s = s*exp(eta*delt);
    for i = 1:N
      for j = 1:3 
        pn(i,j) = pn(i,j)*(1-(eta+epsilon_bar(1,j))*delt/2)/(1+(eta+epsilon_bar(1,j))*delt/2) + F(i,j)*delt/(1+(eta+epsilon_bar(1,j))*delt/2);
      end;
    end;
  end;
  
  % Output the data for every 10 time steps.
  if (tstep >= tstart)
    if (mod(tstep,tramp) == 0)
      position_xyz(tstep, N, partiden, xn, posFile);
      velocity_xyz(tstep, N, partiden, m, pn, velFile);
      totenergy(tstep, N, xn, pn, m, Lo2, eta, Q, Tset, s, a, pressure, estar, kB, rstar, energy_cutoff, totEnergyFile);
      heatcurrent(xn, pn, tstep, Lo2, m, a, partiden, N, energy_cutoff, heatCurrentFile);
    end;
  end;
end;          % End of Verelet algorithm while loop.

% Informing user that program is done.
fprintf('Done!\n');

% Closing all input and output files.
fclose(inFileInitial);
fclose(inFilePhi);
fclose(inFileA);
fclose(posFile);
fclose(velFile);
fclose(totEnergyFile);
fclose(heatCurrentFile);
% ***************************************************************************************************************************************************