units real
  boundary p p p
  atom_style full
 
  
  pair_style lj/cut/coul/long 10.0
  kspace_style pppm 1.0e-4
  bond_style harmonic
  angle_style harmonic 
  dihedral_style opls
  improper_style harmonic
  pair_modify mix geometric tail yes
  special_bonds lj/coul 0.0 0.0 0.5 

  
  
  read_data BTMA_water2.data extra/atom/types 4 extra/bond/types 2 extra/angle/types 2 
  mass 17*18 12.011000 
  mass 19*20 15.999400 
  pair_coeff 17*20 * 0.0 0.0
  bond_coeff 13*14  0.0 0.0
  angle_coeff 17*18 0.0 0.0

  write_data Test.data
 
  

  neighbor 1.5 bin 
  neigh_modify every 5 delay 5 check yes 




velocity all create 50.0 23446733 rot yes mom yes dist gaussian 

   group hydroxide type 7 14
   group interaction type 7 14 1 3 4 5 8 9 12


  
  fix 1 all nve/limit 0.2
  
  run 40000
  
  unfix 1 

  
  minimize 1.0e-6 1.0e-8 1000 10000

  
 
  
  
  fix 2 all nvt temp 50.0 300.0 100 

  thermo_style custom pe density temp press
  thermo 500
  
  run 100000

   unfix 2

   write_restart AfterNVT.restart

   group hydroxide type 7 14

   fix 3 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0

   
   thermo_style custom pe density temp press
   thermo 1000
 
   run 200000

   compute Epair interaction pair/local eng 

   compute Epar all reduce ave c_Epair 
   compute Emin all reduce min c_Epair 

   compute myRDF all rdf 100 4 14 3 14 14 15 
   fix Rdf1 all ave/time 100 100 10000 c_myRDF file tmp.rdf mode vector ave running

   

    fix 4 interaction ave/histo 100 100 10000 -50 180 80 c_Epair file pairE.histo mode vector ave running 

    
    fix 7 all ave/time 100 100 10000 c_Epar  file Epair  
    fix 9 all ave/time 100 100 10000 c_Emin  file Emin  
   

   
   run 2000000

   

   write_restart AfterNPT.restart