#fitting code lammps mohammadreza
# ---------- Initialize Simulation --------------------- 
clear 
units           real 
dimension       3 
boundary        p p p 
  
atom_style      full  
bond_style      class2 
angle_style     class2 
dihedral_style  class2 
improper_style  class2 
# ---------- Define Interatomic Potential --------------------- 
pair_style     soft 10
pair_style lj/class2/coul/long 15.0 

newton        on
comm_modify    mode single cutoff 10 vel no

box tilt large
dielectric 4.5

read_data       hntpu.data

set	atom	1	charge	1.449
set	atom	2	charge	1.449
set	atom	3	charge	1.1
set	atom	4	charge	1.1
set	atom	5	charge	0.2
set	atom	6	charge	0.2
set	atom	7	charge	-0.758
set	atom	8	charge	-0.758
set	atom	9	charge	-0.55
set	atom	10	charge	-0.55
set	atom	11	charge	-0.55
set	atom	12	charge	-0.683
set	atom	13	charge	-0.683
set	atom	14	charge	-0.683
set	atom	15	charge	-0.683
set	atom	16	charge	0.081
set	atom	17	charge	0.625
set	atom	18	charge	-0.727
set	atom	19	charge	-0.272
set	atom	20	charge	-1.2
set	atom	21	charge	-0.127
set	atom	22	charge	0.612
set	atom	23	charge	-0.5
set	atom	24	charge	0.053
set	atom	25	charge	0.353
set	atom	26	charge	-1



##---------------interatomic pot-------------------------------------
kspace_style       pppm 0.00001
pair_modify     mix arithmetic
neighbor        0.5 bin 

special_bonds lj/coul 0.9 0.5 0.75 angle yes dihedral yes

# ---------- define atoms in a small region---------------------
group HNT id < 307
group PU id > 306
 
# ---------- Define Settings -------------------------------------------- 
compute pai all pair lj/class2/coul/long ecoul 
compute        paihnt HNT pair lj/class2/coul/long ecoul 
compute        paipu PU pair lj/class2/coul/long ecoul 
# -------------------------------  
variable y7 equal "c_paihnt"
variable y8 equal "c_paipu"
variable y9 equal "eng"
variable y10 equal (c_pai-(c_paipu+c_paihnt))
# -------------------------------  
variable h7 equal "c_enghnt"
variable h8 equal "c_engpu"
variable h9 equal "eng"
variable h10 equal (c_eng-(c_engpu+c_enghnt))

# ---------- Define Settings ke --------------------------------------------
compute        keratom all ke/atom
compute        ke all reduce sum c_keratom
compute        kehnt HNT reduce sum c_keratom
compute        kepu PU reduce sum c_keratom
# -------------------------------  
variable kk7 equal "c_kehnt"
variable kk8 equal "c_kepu"
variable kk9 equal "ke"
variable kk10 equal (c_ke-(c_kepu+c_kehnt))
# ---------- Define Settings --------------------------------------------
compute eng2 all pe
compute        peratom all pe/atom
compute        pe all reduce sum c_peratom
compute        pehnt HNT reduce sum c_peratom
compute        pepu PU reduce sum c_peratom
# -------------------------------  

variable t7 equal "c_pehnt"
variable t8 equal "c_pepu"
variable t9 equal "pe"
variable t10 equal (c_pe-(c_pepu+c_pehnt))
thermo_style	custom step temp c_pehnt c_pepu etotal press pe c_kehnt c_kepu ke c_pai c_paipu c_paihnt

##---------------ENERGY MINIMIZATION-------------------------------------
fix 15 all box/relax iso 0.0 vmax 0.001

min_style cg
min_modify dmax 0.4
minimize 1e-10 1e-10 10000 100000
reset_timestep 0

run 0 
# ---------- output configuration --------------------- 

variable natoms equal count(all)
variable eperatom equal pe/v_natoms
variable volperatom equal vol/v_natoms
variable teng equal "c_eatoms"
variable t2 equal "epair"
variable t3 equal "ebond"
variable t4 equal "eangle"
variable t5 equal "edihed"
variable t6 equal "etotal"
variable kk equal "ke"
# -------------------------------  
variable h7 equal "c_enghnt"
variable h8 equal "c_engpu"
variable h9 equal "etotal"
#variable h10 equal (c_eng-(c_engpu+c_enghnt))
# -------------------------------  
variable kk7 equal "c_kehnt"
variable kk8 equal "c_kepu"
variable kk9 equal "ke"
#variable kk10 equal (c_ke-(c_kepu+c_kehnt))
# -------------------------------  
variable y7 equal "c_paihnt"
variable y8 equal "c_paipu"
variable y9 equal "c_pai"
variable yy10 equal (c_pai-(c_paipu+c_paihnt))
# ---------- output configuration --------------------- 
print "potential energy per atom = ${eperatom}" 
print "volume per atom = ${volperatom}"

print "pair energy = ${t2}"
print "ebond = ${t3}"
print "eangle = ${t4}"
print "edihed = ${t5}"
print "potential = ${t9}" 
print "potential energy of hnt = ${t7}"
print "potential energy of pu = ${t8}"

print "kinetic energy  = ${kk9}"
print "kinetic energy hnt = ${kk7}"
print "kinetic energy pu = ${kk8}"

variable tpe equal (v_t9-(v_t7+v_t8))
variable tke equal (v_kk9-(v_kk7+v_kk8))

variable int equal (v_tpe+v_tke)

print "total pair energy = ${y9}"
print "pair energy of hnt = ${y7}"
print "pair energy of pu = ${y8}"


#variable h10 equal (v_h9-(v_h7+v_h8))

print "%%ecoh = ${yy10};"
#print "%%ecoh1 = ${h10};"
#print "%%ecoh2 = ${t11};"
print "%%ecoh3 = ${int};"
print "All done"

# output configuration