#!/usr/bin/env python -i
# preceding line should have path for Python on your machine

# Serial syntax: simple.py in.lammps
#                in.lammps = LAMMPS input script

# Parallel syntax: mpirun -np 4 simple.py in.lammps
#                  in.lammps = LAMMPS input script

from __future__ import print_function
import sys
import os
import shutil
import time
import math
import numpy as np
from tabulate import tabulate
#np.set_printoptions(threshold=np.nan)
from scipy import spatial
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

import extfunctions as function

import ctypes

me = 0
from mpi4py import MPI
comm = MPI.COMM_WORLD
me = comm.Get_rank()
nprocs = comm.Get_size()

#Parameters / units
#D_R = 7.82 #microns
#D_M = 50 #vesicle diameter
sigma_r = 0.005 #microns 0.00445 #D_R/D_M #membrane thickness is meant to be about 5nm, where the vesicle diameter is then 50*5nm
sigma_t = 0.01 #microsec #0.025
sigma_m = 1.25E-22 #kg/5nm^3
sigma_v = sigma_r/sigma_t #m/s
sigma_a = sigma_r/(sigma_t**2) #m/s^2
sigma_eta = 1 #sigma_m/(sigma_r*sigma_t) #kg/m.s
d0_M = 1.0 #1.5 #average interparticle separation
d0_R = d0_M*sigma_r
c0_R = -0.562 #microns^-1 										#literature spontaneous curvature of an RBC
c0_M = c0_R*sigma_r
rho_water = 1/(math.pi*4/3)

#thermodynamics
ts = 0.02 #0.005 # 0.05\sigma_t = 0.05x10^-7 sec						# Timestep unit
ini_T = 0.02		                                     	# Sets the initial temperature used in the fix commands below
T_M = 0.23 #0.23 #k_b*T / epsilon	                               	# Sets the final temperature used in the fix commands below
T_damp = 100*ts #1.0                                                   	# Sets the temperature damping parameter used in fix commands below                                                                                       
P_damp = 1.0 #1000*ts                                             	# Sets the Pressure damping variable used in the "fix npt" commands
P = 0.05                                                  	# Sets the desired system pressure in the "fix npt" commands

epsilon_lj = 0.2                                                   	# Sets the energy well (epsilon) variable used in the pair coefficient specifications
sigma_lj = 2.7

#Pair-potential parameters 	
rcut_yuan = 2.6	#2.7 #2.8 #2.6												# Pair-potential cutoff length
epsilon_yuan = 1.5 #1.5 #1.1
sigma_yuan = 1.0
zeta_yuan = 4 #4
mu_yuan = 3 #6 #3
beta_R = (c0_R*d0_R)/2 #0 #0.0194*0.0		#beta = sin(theta_0) =c_0*d_0/2
beta_M = (c0_M*d0_M)/2 #0


def gather_atoms():
	gath_pos = lmp.gather_atoms_concat("x", 1, 3) # (name,type,count) return per-atom property of all atoms gathered into data, ordered by atom ID
                                          # name = "x", "charge", "type", etc
	gath_pos = np.array(gath_pos).reshape(-1,3)
	gath_type = lmp.gather_atoms_concat("type", 0, 1)
	gath_id = lmp.gather_atoms_concat("id", 0, 1)
	gath = np.column_stack((gath_id, gath_type, gath_pos))

	#print(gath.shape)

	return gath

def get_nwater():
	lmp.command('variable nwater equal count(water_in)')
	#comm.Barrier()
	n_water = lmp.extract_variable('nwater', 'all', 0)  # extract value(s) from a variable
	#comm.Barrier()
	nwdel = n_water - ramp2*n_water
	#comm.Barrier()
	if me == 0:
		print("Water_in = {}, nwdel = {}".format(n_water, nwdel))
	
	return nwdel

def main():
	#lmp.command('package intel 0 omp 2 mode double')
	#lmp.command('package intel 1')
	lmp.command('suffix opt') #suffix off
	#lmp.command('suffix hybrid intel opt')
	
	# ------------------------ VARIABLES ----------------------------
	lmp.command('variable 		scale1 equal ramp({},{})'.format(ramp[0], ramp[1])) #2.7, 2.52  # Sets the Lj potential's sigma value in the final water fix to decay from 2.7 to 2.52
	lmp.command('variable 		logfid string log/log.{}'.format(fid))


	# ------------------------ INITIALIZATION ----------------------------
	lmp.command('units			lj')                            			# Sets the units used for the simulation, in this case: Lennard-Jones units                                                       
                			
	lmp.command('atom_style    	hybrid 		ellipsoid peri molecular')		# Defines the style used for the atoms, where each style allows atoms to have certain properties. 
								#"hybrid" allows combination of multiple styles, to give access to each of there associated paramters 
											#ellipsoid = for aspherical particles; (shape, quaternion, angular momentum)
											#peri = for spherical, mesocopic Peridynamic models; (mass, volume). Used here for extracting the initial configurations (x0) of all the CG particles
											#molecular = for uncharged molecules; (bonds, angles, dihedrals, impropers)
	#lmp.command('atom_style   	hybrid 		ellipsoid molecular')

	lmp.command('bond_style   	harmonic1')                     			# Defineds how bonds are calculated, where bonded atoms are defined as such in the data file
								#harmonic1 is a custom style, for the harmonic spring defined by Hooke's law. This is used for the bonding between the bilayer membrane and cytoskeleton

	lmp.command('pair_style   	hybrid   lj/cut 3.6  fluidmembrane {}'.format(rcut_yuan))    # Sets the pair_style potentials used for the simulation, in the case both lennard-jones/cut and the fluidmembrane potential
										#L-J is used for interactions between CG lipid, cytoskeleton and water particles 
	lmp.command('boundary     p  p  p')                                     # Creates the periodic boundaries of the simulation box

	box_info = lmp.extract_box()
	if me == 0:
		print('Prediodicity = ', box_info[5])                                                

	#comm.Barrier()
	#print("Proc %d out of %d procs has" % (me,nprocs), lmp)
	#comm.Barrier()

	# ------------------------ ATOM DEFINITION/SETTINGS ----------------------------
	lmp.command('read_data    	{}'.format(datafile))                      					# Specifies which read data file will be used to read in initial coordinates of system
	proc_dir = 'proc/proc.{}_runs{}{}'.format(fid, run_T, note)
	#lmp.command("processors * * * file outfile {}".format(proc_dir))	

	#lmp.command('mass 		*	1.0')
	lmp.command('set type       	* mass     		1.0')
	#lmp.command('set type       	1*5 mass     		1.0')      										# Sets mass of all particles to 1.0 (should it be just the bilayer that's set?)                                           
	#lmp.command('set type       	6 	mass     		{}'.format(4*rho_water))      										# Sets mass of all particles to 1.0 (should it be just the bilayer that's set?)                                           
	#lmp.command('set type       	7 	mass 			{}'.format(4*rho_water))      										# Sets mass of all particles to 1.0 (should it be just the bilayer that's set?)                                           
	#lmp.command('set type       	8   mass 			1.0')      										# Sets mass of all particles to 1.0 (should it be just the bilayer that's set?)                                           

	lmp.command('set type     	1*2 shape  	1 1 1')   					     				# membrane particle are rigid body (ellipsoids), needing shape to be defined

	lmp.command('group        	bilayer  	type 1 2')                                      # Creates a group called bilayer which consists of type 1 and 2 particles, (Type 1 Bilayer, Type 2 Transmembrane proteins)
	lmp.command('group        	network   	type 3 4 5')                                    # Creates a group called network which consists of type 3 4 5 particles
	lmp.command('group        anchor_bonds  type 2 3 5')
	lmp.command('group        	water_in  	type 6')                                        # Creates a group called water_in consisting of type 6 particles ( Type 6 = Water inside Membrane)
	lmp.command('group        	water_out 	type 7')                                        # Creates a group called water_out consisting of type 7 particles ( Type 7 = Water Outside Membrane)
	lmp.command('group        	water   	type 6 7')                                      # Creates a group called water consisting of both type 6 and 7 particles ( Water inside + Out)
	lmp.command('group        	tube   		type 8')                                      # Creates a group called tube consisting of the tube boundary particles
	#lmp.command('group        	cell   		type 1 2 3 4 5 8')                                      # Creates a group called tube consisting of the tube boundary particles
	#lmp.command('group        	cell   		type 1 2 3 4 5 6 7 8')                                      # Creates a group called tube consisting of the tube boundary particles
	lmp.command('group        	cell   		type 1 2')                                      # Creates a group called tube consisting of the tube boundary particles

	lmp.command('velocity     	bilayer 	create {} 87287 loop geom'.format(T_M))           # Creates Initial Velocity for bilayer particles ( Type 1 and 2)
	lmp.command('velocity     	network 	create {} 87287 loop geom'.format(ini_T))       # Creates Initial Velocity for network particles ( Type 3 and 4 and 5)
	lmp.command('velocity     	water 		create {} 87287 loop geom'.format(T_M))           # Creates Initial Velocity for water particles ( Type 6 and 7 )
	#lmp.command('velocity     	tube 		zero linear')           # Creates Initial Velocity for water particles ( Type 6 and 7 )
	#lmp.command('region 		boundary 	cylinder x 30 30 60 -60.0 120.0')

	lmp.command('neighbor    	1.0 		bin') #Creates bin-list of all atoms neighbouring eachother by 1.0
	
	lmp.command('region rsphere block INF INF INF INF INF INF') #sphere 0 0 0 100') #Creates bin-list of all atoms neighbouring eachother by 1.0
	          
	# ------------------------ FORCE FIELDS ------------------------------	
	
	lmp.command('bond_coeff   	1			50  	1.5')         	# Defines the bond coefficients for the spectrin network bonds 
								#bond-type, K, 		r0
	lmp.command('bond_coeff   	2*3 		20		2.0')  			# Defines the bond coefficients for the membrane-spectrin anchoring 
								#bond-type, K, 		r0

	lmp.command('pair_coeff   	1*2 	1*2 	fluidmembrane  	{}  {}  {}  {}  {}  {}'.format(epsilon_yuan, sigma_yuan, rcut_yuan, zeta_yuan, mu_yuan, beta_M))   # Sets the pair coefficients and parameters between type 1 and 2 (bilayer) particles
	lmp.command('pair_coeff   	1*2    	3*5   	lj/cut   		{} 		{}'.format(epsilon_lj, sigma_yuan))  	# Sets L-J equilibirium length (sigma) and constant epsilon for type 1,2 (bilayer) and type 3,4,5 (network) particle interactions 
	lmp.command('pair_coeff   	1*2    	6   	lj/cut   		{} 		{}'.format(epsilonWI_lj, sigmaW_yuan))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 1,2 (bilayer) and type 6,7 (water) particle interactions
	lmp.command('pair_coeff   	1*2    	7   	lj/cut   		{} 		{}'.format(epsilonWE_lj, sigmaW_yuan))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 1,2 (bilayer) and type 6,7 (water) particle interactions
	lmp.command('pair_coeff   	1*2    	8   	lj/cut   		{} 		{}'.format(epsilon_lj, sigma_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 1,2 (bilayer) and type 8 (tube) particle interactions
	
	lmp.command('pair_coeff   	3*5    	3*5   	lj/cut   		{} 		{}'.format(epsilon_lj, sigma_yuan))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 3,4,5 (network) and type 3,4,5 (network) particle interactions
	lmp.command('pair_coeff   	3*5    	6   	lj/cut   		{} 		{}'.format(epsilonWI_lj, sigmaW_yuan))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 3,4,5 (network) and type 6,7 (water) particle interactions
	lmp.command('pair_coeff   	3*5    	7   	lj/cut   		{} 		{}'.format(epsilonWE_lj, sigmaW_yuan))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 3,4,5 (network) and type 6,7 (water) particle interactions
	lmp.command('pair_coeff   	3*5    	8   	lj/cut   		{} 		{}'.format(epsilon_lj, sigma_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 8 (boundary) and type 8 (tube) particle interactions                                                       
	
	lmp.command('pair_coeff   	6    	7   	lj/cut   		{} 		{}'.format(epsilonWE_lj, sigmaW_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 6,7 (water) and type 6,7 (water) particle interactions                                                       
	lmp.command('pair_coeff   	6    	6   	lj/cut   		{} 		{}'.format(epsilonWI_lj, sigmaW_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 6,7 (water) and type 6,7 (water) particle interactions                                                       
	lmp.command('pair_coeff   	7    	7   	lj/cut   		{} 		{}'.format(epsilonWE_lj, sigmaW_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 6,7 (water) and type 6,7 (water) particle interactions                                                       
	lmp.command('pair_coeff   	6*7    	8   	lj/cut   		{} 		{}'.format(epsilonWE_lj, sigmaW_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 6,7 (water) and type 8 (tube) particle interactions                                                       
	
	lmp.command('pair_coeff   	8    	8   	lj/cut   		{} 		{}'.format(epsilon_lj, sigma_lj))    # Sets L-J equilibirium length (sigma) and constant epsilon for type 8 (tube) and type 8 (tube) particle interactions                                                       
	#sigma for water-water interactions controls the internal RBC volume (smaller sigma means smaller effective volume for the water)
	

	# ------------------------ OUTPUT ------------------------------
	
	lmp.command("variable mass equal mass(all)/count(all)")
	lmp.command("variable mass_bilayer equal mass(bilayer)/count(bilayer)")
	lmp.command("variable mass_network equal mass(network)/count(network)")
	lmp.command("variable mass_waterin equal mass(water_in)/count(water_in)")
	lmp.command("variable mass_waterout equal mass(water_out)/count(water_out)")
	lmp.command('thermo_style 	custom step temp press v_mass')   	# Printing thermodynamic data to the screen and log file 

	lmp.command('thermo	     	200')                               # Sets the number of timesteps to run before printing out thermodynamic data 
	lmp.command('timestep     	{}'.format(ts))					# Sets the timestep \Delta.t of the simulation
	lmp.command('thermo_modify lost warn')

	#lmp.command('compute sep bilayer coord/atom cutoff {} 1*2'.format(rcut_yuan))

	#lmp.command('variable X atom x')
	lmp.command('variable r atom sqrt(y*y+z*z)')
	lmp.command('variable FX atom v_F*2*(1-((v_r*v_r)/({}*{})))'.format(R_tube, R_tube)) #4*step

	lmp.command("variable    p equal 400")     # correlation length
	lmp.command("variable    s equal 5")       # sample interval
	lmp.command("variable    d equal $p*$s")   # dump interval 
	lmp.command("variable scale1 equal (1.0/{})*vol*$s*dt".format(T_M))
	lmp.command("variable scale2 equal vol/(2.0*{}*dt*$d)".format(T_M))
	lmp.command("variable pxy equal pxy")
	lmp.command("variable pxz equal pxz")
	lmp.command("variable pyz equal pyz")

	dump_dir = 'dump/dump.{}_runs{}{}.lammpstrj'.format(fid, run_T, note)    
	write_dir = 'data/data.{}_runs{}{}'.format(fid, run_T, note)
	writer_dir = 'data/restart.{}_runs{}{}'.format(fid, run_T, note)
	
	#lmp.command('dump 1 all atom 2000 {}'.format(dump_dir))          # Command sets when to dump a snapshot of atom quantities to one or more files every N timesteps in one of several styles.
	lmp.command('dump 1 cell atom 2000 {}'.format(dump_dir))          # Command sets when to dump a snapshot of atom quantities to one or more files every N timesteps in one of several styles.

	#mov_dir = 'dump/mov.{}_runs{}{}.mpg'.format(fid, run_T, note)                                                          
	#lmp.command('dump 2 cell movie 1000 {} type type size 1280 720'.format(mov_dir))
		
	#dump the final atom positions of the bilayer and cytoskeleton	
	#lmp.command('dump report_all all custom {} all.report id type mass'.format(run_T))
	#lmp.command('dump_modify report_all delay {}'.format(run_equilib))
	# lmp.command('dump report_bilayer bilayer custom {} bilayer.report id type x y z'.format(run_T))
	# lmp.command('dump_modify report_bilayer delay {}'.format(run_equilib))
	# lmp.command('dump report_network network custom {} network.report id type x y z'.format(run_T))
	# lmp.command('dump_modify report_network delay {}'.format(run_equilib))
	# lmp.command('dump report_forces all custom {} forces.report id type fx fy fz'.format(run_T))	
	# lmp.command('dump_modify report_forces delay {}'.format(run_equilib))
	
	#lmp.command('dump report_mass all custom {} mass.report id type mass'.format(run_T))
	#lmp.command('dump_modify report_mass delay {}'.format(run0))

	#lmp.command("balance 1.0 x uniform y uniform z uniform out {}".format(proc_dir)) ##default
	#lmp.command("balance 1.0 shift x 500 1.0 out {}".format(proc_dir))
	lmp.command("fix balance all balance 5000 1.0 shift xyz 20 1.0 weight group 2 water 0.5 bilayer 6.0 out {}".format(proc_dir))

	# ------------------------ SIMULATION ----------------------------- #
	start = time.time()
	timerR = 0
	natoms = lmp.get_natoms()

	#gath = gather_atoms()
	comm.Barrier()
	headers = ["Elapsed[s]", "Step", "Time[micro-s]", "V[microns^3]", "A[microns^2]", "D_major[r^M]", "D_minor[r^M]", "V/D", "viscosity"]
	if me == 0:
		print('\nWILL BE DUMPED TO: ', dump_dir)
		VA = np.empty([0, 9])

	burst = False
	
	mass_bilayer = lmp.extract_variable('mass_bilayer', 'all', 0)
	mass_network = lmp.extract_variable('mass_network', 'all', 0)
	mass_waterin = lmp.extract_variable('mass_waterin', 'all', 0)
	mass_waterout = lmp.extract_variable('mass_waterout', 'all', 0)
	
	if me == 0:
		log_dir = 'log/log_{}_runs{}{}.txt'.format(fid, run_T, note)
		if os.path.isfile(log_dir):
			os.remove(log_dir)
			print("DELETED OLD log FILE ", log_dir)
		f = open(log_dir, 'w+')
		f.write('T_M = {}, sigma_r = {}, d0_M = {}, ts = {}, sigma_t = {}, c0_R = {}\n'.format(T_M, sigma_r, d0_M, ts, sigma_t, c0_R))
		f.write('P_damp = {}, P = {}, ini_T = {}, T_damp = {}\n'.format(P_damp, P, ini_T, T_damp))
		f.write('ramp = ({},{}), ramp2 = {}, epsilon_lj = {}, sigma_lj = {}, epsilonWI_lj = {}, epsilonWE_lj = {}, sigmaW_lj = {}, force = [{},{},{}]\n'.format(ramp[0], ramp[1], ramp2, epsilon_lj, sigma_lj, epsilonWI_lj, epsilonWE_lj, sigmaW_lj, *F))
		f.write('rcut_yuan = {}, epsilon_yuan = {}, sigma_yuan = {}, zeta_yuan = {}, mu_yuan = {}, beta_M = {:.5f}, beta_R = {:.5f}\n'.format(rcut_yuan, epsilon_yuan, sigma_yuan, zeta_yuan, mu_yuan, beta_M, beta_R))		
		f.write('D_M = {}\n'.format(D_M))
		f.write('mass_bilayer = {:.4f}, mass_network = {:.4f}, mass_waterin = {:.4f}, mass_waterout = {:.4f}\n'.format(mass_bilayer, mass_network, mass_waterin, mass_waterout))


		f.write('\n{}, {}, {}, {}, {}, {}, {}, {}, {}\n'.format(*headers))
		f.close()

	for t in range(0, run_T, run_dt):
		comm.Barrier()
		if (burst == True):
			break
		t1 = time.time()

#Tube equilibriate
		if (t < run0_T):
			if (t == 0):
				c0 = 'fix 0 tube freeze'
				lmp.command(c0)
			if (me==0): 
				print('\n', c0)
				print('run {}/{}\n'.format(t+run_dt, dt_0*run_dt))		#int(run_dt/10)))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, 0, run0_T))		#int(run_dt/10)))	#10000
			#lmp.command('info all out overwrite log/log_info.txt')	#10000

#Water equilibriate
		elif (t>=run0_T and t<run1_T):
			if (t == run0_T):
				c1 = 'fix w1 water_out npt temp {} {} {} iso {} {} {}'.format(ini_T, T_M, T_damp, P, P, P_damp)
				#c1_2 = 'fix w2 water_in npt temp {} {} {} iso {} {} {}'.format(T_M, T_M, T_damp, P, P, P_damp)
				#c1 = 'fix w1 water_out nvt temp {} {} {}'.format(ini_T, T_M, T_damp)
				c1_2 = 'fix w2 water_in nvt temp {} {} {}'.format(ini_T, T_M, T_damp)
				lmp.command(c1) # Applies an NPT ensemble fix to the water group 
				lmp.command(c1_2) # Applies an NPT ensemble fix to the water group 

			if (me==0): 
				print('\n', c1)
				print('\n', c1_2)
				print('run {}/{}\n'.format(t+run_dt-run0_T, dt_1*run_dt))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, run0_T, run1_T))	#10000
			#lmp.command('info all out append log/log_info.txt')	#10000

#Network equilibriate
		elif (t>=run1_T and t<run2_T):
			if (t == run1_T):
				c2 = 'fix 2 network npt temp {} {} {} iso {} {} {}'.format(ini_T, T_M, T_damp, P, P, P_damp)
				lmp.command(c2) # Applies an NPT ensemble fix the network group using same parameters as above
				#perform time integration on Nose-Hoover style non-Hamiltonian equations of motion, calculating positions and velocities sampled from the isothermal-isobaric (npt) ensemble.
			if (me==0): 
				print('\n', c2)
				print('run {}/{}\n'.format(t+run_dt-run1_T, dt_2*run_dt))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, run1_T, run2_T))	#20000
			#lmp.command('info all out append log/log_info.txt')	#10000

#Bilayer equilibriate
		elif (t>=run2_T and t<run3_T):
			if (t == run2_T):
				c3 = 'fix 3 bilayer nvt/asphere temp {} {} {}'.format(ini_T, T_M, T_damp)
				lmp.command(c3) # Applies an NVT (not NPT) ensemble fix to the network group using only temperature parameters               
			#"fix nvt/asphere" performs constant NVT integration to update position, velocity, orientation, and angular velocity each timestep for aspherical or ellipsoidal particles in the group using a Nose/Hoover temperature thermostat. 
			#"fix nvt" command instead assumes point particles and only updates their position and velocity.
			if (me==0): 
				print('\n', c3)
				print('run {}/{}\n'.format(t+run_dt-run2_T, dt_3*run_dt))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, run2_T, run3_T)) 	#10000
			#lmp.command('info all out append log/log_info.txt')	#10000

#Water condense
		elif (t>=run3_T and t<run4_T):
			if (t == run3_T):
				gath = gather_atoms()

				squish_T, squish_B, n_bilayer = function.squish(gath, sigma_r, R_squish)
				
				# if me == 0:
					# print("n_squishT = {} ({:.2f}%), n_squishB = {} ({:.2f}%)".format(len(squish_T), 100*(len(squish_T)/n_bilayer), len(squish_B), 100*(len(squish_B)/n_bilayer)))
				# squish_T = squish_T[:,0].astype(int)
				# str_T = ' '.join(str(e) for e in squish_T)
				# squish_B = squish_B[:,0].astype(int)
				# str_B = ' '.join(str(e) for e in squish_B)
				# lmp.command("group squishT id {}".format(str_T))
				# lmp.command("group squishB id {}".format(str_B))

				lmp.command('region	squishT	cylinder z 0 0 {} 0 INF'.format(R_squish))
				lmp.command('region	squishB	cylinder z 0 0 {} INF 0'.format(R_squish))

				#c4 = 'fix 4 water_in adapt 1 pair lj/cut sigma 6 6 v_scale1'
				nwdel = get_nwater()
				c4 = 'fix 4 water_in evaporate {} {} rsphere 38277 molecule no'.format(run_dt, int(round(nwdel/dt_4))) #N = 1000 delete atoms every this many timesteps M = 10 number of atoms to delete each time
				lmp.command(c4) 	# Applies an adaptive (change over-time) fix to the LJ potential's sigma for the water group, having it decay over time according to scale1 

			frac = [1/len(squish_T),1/len(squish_B)]
			#c5 = 'fix 5 squishT addforce 0 0 {} every 1'.format(frac[0]*scale*F_squish[0]) #-20
			#c6 = 'fix 6 squishB addforce 0 0 {} every 1'.format(frac[1]*scale*F_squish[1]) #20
			c5 = 'fix 5 bilayer addforce 0 0 {} every 1 region squishT'.format(frac[0]*scale*F_squish[0]) #-20
			c6 = 'fix 6 bilayer addforce 0 0 {} every 1 region squishB'.format(frac[1]*scale*F_squish[1]) #20
			lmp.command(c5)
			lmp.command(c6)

			if (me==0): 
				print(c4)
				print(c5)
				print(c6)
				print('run {}/{}\n'.format(t+run_dt-run3_T, dt_4*run_dt))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, run3_T, run4_T))	#40000									# This effectively reduces the internal water volume
			#lmp.command('info all out append log/log_info.txt')	#10000

#Simulate flow
		elif (t>=run4_T and t<run5_T):
			if (t == run4_T):
				gath = gather_atoms()
				if me == 0:
					D0_major, D0_minor, D0_z = function.get_Dsimp(gath, sigma_r)
					print('\nWidth of RBC = {:.2f}[r^M], {:.2}[microns]'.format(D0_major/sigma_r, D0_major))
				
				lmp.command('unfix 4')
				lmp.command('unfix 5')
				lmp.command('unfix 6')
				lmp.command("fix SS all 		ave/correlate $s 		$p 		$d 		v_pxy v_pxz v_pyz type auto ave running") #file S0St.dat ave running")
							#fix ID group-ID 	ave/correlate Nevery 	Nrepeat Nfreq 	value1 value2 ... keyword args ...
				lmp.command("variable v11 equal trap(f_SS[3])*${scale1}")
				lmp.command("variable v22 equal trap(f_SS[4])*${scale1}")
				lmp.command("variable v33 equal trap(f_SS[5])*${scale1}")
				#lmp.command("variable vgk equal (v_v11+v_v22+v_v33)/3.0")
				
				c5 = 'fix 5 water_out addforce v_F 0 0 every {}'.format(1) #v_F #v_FX #*F
				lmp.command(c5) #run_dt/100
			if (me==0): 
				print('\nunfix 4')
				print('\n', c5)
				print('run {}/{}\n'.format(t+run_dt-run4_T, dt_5*run_dt))
			lmp.command('run {} start {} stop {} post no'.format(run_dt, run4_T, run5_T))	#40000									# This effectively reduces the internal water volume
			#lmp.command('info all out append log/log_info.txt')	#10000
		
		comm.Barrier()
		t2 = time.time()
		timerR += (t2 - t1)

		#gath = gather_atoms()
		comm.Barrier()

		nwdel = get_nwater()

		if (t > run4_T):
			#lmp.command("variable eta_ein equal {}*(v_vxy+v_vxz+v_vyz)/3.0".format(sigma_eta))
			#lmp.command("print 'average viscosity (Einstein): ${eta_ein}'")
			#viscEIN_M = lmp.extract_variable('eta_ein', 'all', 0)
			
			lmp.command("variable eta_gk equal {}*(v_v11+v_v22+v_v33)/3.0".format(sigma_eta))
			lmp.command('print "average viscosity (GK): ${eta_gk}"')
			visc_M = lmp.extract_variable('eta_gk', 'all', 0)
		else:
			#viscEIN_M = 0.0
			visc_M = 0.0

		if (me == 0): #and (t+run_dt)%(2*run_dt)==0):
		
			step = lmp.get_thermo('step')
			t3 = time.time()
			V, A = 0.0, 0.0 #function.get_VA(gath, sigma_r)		
			t4 = time.time()
			D_major, D_minor, D_z = 0.0, 0.0, 0.0 #function.get_Dsimp(gath, sigma_r)
			t5 = time.time()
			R = (D_major + D_minor)/4
			
			if (R != 0): 
				V_coef = V/(R**3)
				V0 = (2/3)*math.pi*R**3
				V_coef2 = V/V0
			else: 
				V_coef = 0.0
				V_coef2 = 0.0
			A_coef = 0.0 # A/V
			print("V_coef = {:.2f} ({:.1f}%, 1.57), V_coef2 = {:.2f} ({:.1f}%, 0.64), A_coef = {:.2f} ({:.1f}%, 1.44)".format(V_coef, 100*(V_coef/1.57 - 1), V_coef2, 100*(V_coef2/0.64 - 1), A_coef, 100*(A_coef/1.44 - 1)))
				
			VA = np.vstack([VA, [t2-t1, int(step), step*sigma_t, V, A, D_major/sigma_r, D_minor/sigma_r, V_coef, visc_M]])
			print('\n D_major = {:.2f}[r^M], {:.5}[microns]'.format(D_major/sigma_r, D_major))
			print("Average viscosity = {}".format(visc_M))
			print('(run took {:.3f}[s]) Step {:.0f} at {:.2f}[micro-sec]: V = {:.5f}[microns^3], A = {:.5f}[microns^2]\n'.format(*VA[len(VA)-1]))
			t_tot = time.time() - t1
			print('TIMINGS: run={:.4f}, get_VA={:.4f}, get_Dsimp={:.4f}'.format((100*(t2-t1)/t_tot), (100*(t4-t3)/t_tot), (100*(t5-t4)/t_tot)))
			print('\n ---------- ELAPSED TIME: {:.3f} [h] ----------'.format((time.time() - start)/3600))
			f = open(log_dir, 'a')
			f.write('{:.3f} 	    {:.0f} 	{} 	{:.5f} 	{:.5f}	{:.5f}	{:.5f}	{:.5f}	{}\n'.format(*VA[len(VA)-1]))
			f.close()

			if ((len(VA)-5 >= 0) and VA[len(VA)-1, 3] > 2*VA[len(VA)-4, 3]):
				print("\n    MEMBRANE BURST at step {}\n".format(t))
				burst = True
				#t = run_T
				f = open(log_dir, 'a')
				f.write("\n    MEMBRANE BURST at step {}\n".format(t))
				f.close()					
			else:
				burst = False
		else:
			burst = False
		burst = comm.bcast(burst, root=0)


	# -----------------------------------------------------------------

	end = time.time()
	
	# lmp.command('unfix 0')
	# lmp.command('unfix w1')
	# lmp.command('unfix w2')
	# lmp.command('unfix 2')
	# lmp.command('unfix 3')
	# lmp.command("write_data {} nocoeff nofix".format(write_dir))
	# lmp.command("write_restart {}".format(writer_dir))
	# if me == 0: 
		# print("WROTE DATA FILE TO {}".format(write_dir))
		# print("WROTE RESTART FILE TO {}".format(writer_dir))

	if me == 0:

		# tabulate data
		table = tabulate(VA, headers, tablefmt="simple")
		print(table)

		print('\n EXECUTION TIME: {:.3f} [h]'.format((end - start)/3600))
		print('\n RUN TIME FRAC: {:.3f} [h]'.format(100*(timerR/(end-start))))
		print(' DUMPED TO: ', dump_dir)
	
		f = open(log_dir, 'a')	
		f.write('\n EXECUTION TIME: %.3f [h]' % ((end - start)/3600))
		f.write('\n RUN TIME FRAC: %.3f [h]' % (100*(timerR/(end-start))))
		f.close()

if __name__ == '__main__':

	argv = sys.argv
	if len(argv) < 4:
		note = ''
		if me == 0:
			print("Syntax: python RBC1_tube_nwdelNN.py 1111111 data.BLAH NOTE")
	elif len(argv) < 3:  	
		sys.exit()
		if me == 0:
			print("Syntax: python RBC1_tube_nwdelNN.py 1111111 data.BLAH NOTE")
	else:
		note = '_' + sys.argv[3]

	fn_py1 = sys.argv[0]
	slurm_id = sys.argv[1]
	slurm_name = 'slurm-' + slurm_id + '.out'
	if me == 0:
		print("SLURM_NAME = ", slurm_name)

	datafile = sys.argv[2]
	D_M = int(datafile[datafile.find("D")+1:datafile.find("A")])
	if (D_M == 62 or D_M == 87 or D_M == 112):
		D_M = D_M + 0.5
	R_tube = int(datafile[datafile.find("@")+1:datafile.find("g")]) + (D_M/2)
	factor_water = int(datafile[datafile.find("W")+1:datafile.find("_@")])
	scale = (D_M/50)**2
	epsilonWI_lj = epsilon_lj #*factor_water
	epsilonWE_lj = epsilon_lj #*factor_water
	factor_water2 = factor_water
	sigmaW_yuan = sigma_yuan*factor_water2
	sigmaW_lj = sigma_lj*factor_water2
	
	ramp2 = 1.57-(D_M/50)**(1/5)
	if (D_M == 50):
		ramp = [sigma_lj*factor_water2, 2.45*factor_water2] #2.55	#2.58 #[2.8, 2.62] #sigma_ew
		ramp2 = 0.68 #0.50 #0.55 #0.38 #0.3
		F_squish = [-0.5, 0.5] #3
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 62):
		ramp = [sigma_lj*factor_water2, 2.45*factor_water2] #2.55	#2.58 #[2.8, 2.62] #sigma_ew
		ramp2 = 0.58 #0.50 #0.55 #0.38 #0.3
		F_squish = [-0.5, 0.5] #6
		R_squish = (0.6)*(D_M/2) #0.5 #D_M/4 #radius of circle to push
	elif (D_M == 75):
		ramp = [sigma_lj*factor_water2, 2.45*factor_water2] #2.55	#2.58 #[2.8, 2.62] #sigma_ew
		ramp2 = 0.54 #0.50 #0.55 #0.38 #0.3
		F_squish = [-1, 1] #6
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 87):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		ramp2 = 0.42 #0.46 #0.47 #0.35 #0.3
		F_squish = [-5, 5] #7
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 100):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.42 #0.46 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 112):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.38 #0.43 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 125):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.38 #0.43 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 150):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.33 #0.41 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.5)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 175):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.33 #0.41 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.5)*(D_M/2) #D_M/4 #radius of circle to push
	elif (D_M == 200):
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.33 #0.41 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.5)*(D_M/2) #D_M/4 #radius of circle to push
	else:
		ramp = [sigma_lj*factor_water2, 2.00*factor_water2] #2.10	#2.58 #[2.8, 2.62] #sigma_ew
		#ramp2 = 0.30 #0.33 #0.41 #0.47 #0.35 #0.3
		F_squish = [-10, 10]
		R_squish = (0.6)*(D_M/2) #D_M/4 #radius of circle to push	
	#F_squish = np.array(F_squish)
	F_squish = [0.0, 0.0]

	if me == 0:
		print("DATAFILE = ", datafile)
		print("R_tube = ", R_tube)
		print("D_M = ", D_M)
		print("factor_water = ", factor_water)
		print("Scale = ", scale)
		print("ramp = [{},{}]".format(ramp[0], ramp[1]))
		print("ramp2 = {}".format(ramp2))
		print("F_squish = ", F_squish)

	fid = datafile[5:]
	if me == 0:
		print("FID = ", fid)

	if me == 0:
		print("\nBeta^R = {}, theta^R = {}".format(beta_R, np.arcsin(beta_R)))
		print("Beta^M = {}, theta^M = {}\n".format(beta_M, np.arcsin(beta_M)))
		
	run_dt = 5000
	dt_0 = 1 				#fix boundary tube
	dt_1 = 5 #2 #1 				#fix water
	dt_2 = 5 #4 #2 #2 			#fix network
	dt_3 = 10 #7 #8 #4 #1			#fix bilayer
	#if (D_M >= 50 and D_M < 75):
	#	dt_4 = 12 #8
	#elif (D_M >= 75 and D_M < 100):
	#	dt_4 = 12 #24 #16 #8 #3			#reduce internal water volume
	#else:
	#	dt_4 = 24
	dt_4 = 12
	dt_5 = 25 #20 #10 #8 #6 	#fix flow

	run0 = dt_0*run_dt; run0_T = run0
	run1 = dt_1*run_dt; run1_T = run0 + run1
	run2 = dt_2*run_dt; run2_T = run0 + run1 + run2
	run3 = dt_3*run_dt; run3_T = run0 + run1 + run2 + run3
	run4 = dt_4*run_dt; run4_T = run0 + run1 + run2 + run3 + run4
	run5 = dt_5*run_dt; run5_T = run0 + run1 + run2 + run3 + run4 + run5

	run_equilib = run0 + run1 + run2 + run3 + run4
	run_T = run_equilib + run5


	from lammps import lammps #, PyLammps
	lmp = lammps()
	#L = PyLammps(ptr=lmp)

	F = [0.0, 0.0, 0.0] #0.04
	#F = [100*0.040, 0, 0] #[0.02, 0, 0] #Force applied to external water particles to create flow
	
	lmp.command('variable F equal {}'.format(F[0])) #4*step
	#F_L = [-XF, -0.0, -0.0] #-0.1 #0.04
	#F_R = [XF, 0.0, 0.0] #0.1 #0.04
	comm.Barrier()
	main()
	#comm.Barrier()
	
	if me == 0:
		if os.path.isfile(slurm_name):
			fn = 'slurm.{}_runs{}{}'.format(fid, run_T, note)
			os.rename(slurm_name, fn)
			print("FOUND ", slurm_name)
		else:
			print("NOT FOUND ", slurm_name)
			
		fn_py2 = 'pyfiles/{}_runs{}{}.py'.format(fid, run_T, note)
		shutil.copy2(fn_py1, fn_py2)

	#MPI.Finalize()