# -*- python script-*-
# MD code of Stinger-Weber Silicon

# Definition of procedures

import mdpp

def initmd(T) :
    mdpp.cmd('''
    setnolog
    setoverwrite
    ''')
    mdpp.cmd('dirname = runs/si-example-{0}'.format(T))

def makecrystal(nx,ny,nz) :
    # create perfect lattice configuration
    mdpp.cmd('crystalstructure = diamond-cubic')
    mdpp.cmd('latticeconst = 5.4309529817532409') # (A) for Si
    mdpp.cmd('element0 = Silicon')
    mdpp.cmd('''
    latticesize = [ 1 0 0 {0}
                    0 1 0 {1}
                    0 0 1 {2} ]
    '''.format(nx,ny,nz))
    mdpp.cmd('makecrystal  finalcnfile = perf.cn  writecn')


def setup_window() :
    # plot settings
    mdpp.cmd('''
     atomradius = 0.67 bondradius = 0.3 bondlength = 2.8285 #for Si 
     atomcolor = orange highlightcolor = purple 
     bondcolor = red backgroundcolor = gray70
    #plot_color_bar = [ 1 -4.85 -4.50 ]  highlightcolor = red
     plotfreq = 10  rotateangles = [ 0 0 0 1.25 ]
    ''')

def openwindow() :
    #Configure and open a plot
    setup_window()

    #### setup_window
    mdpp.cmd('openwin  alloccolors rotate saverot eval plot')
    

def relax_fixbox() :
    #Conjugate-Gradient relaxation

    mdpp.cmd('''
    conj_ftol = 1e-7 conj_itmax = 1000 conj_fevalmax = 10000 
    conj_fixbox = 1 #conj_monitor = 1 conj_summary = 1 
    relax
    ''')
    

def setup_md() :
    # MD settings (without running MD)

    mdpp.cmd('''
    equilsteps = 0  totalsteps = 100 timestep = 0.0001 # (ps) 
    atommass = 28.0855 # (g/mol) 
    DOUBLE_T = 0 
    randseed = 12345 srand48
    initvelocity totalsteps = 10000 saveprop = 0 
    saveprop = 1 savepropfreq = 10 openpropfile 
    vt2=1e28
    ensemble_type = NVE integrator_type = VVerlet 
    implementation_type = 0
    
    totalsteps = 5000 
    output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESS_xx TSTRESS_yy TSTRESS_zz" 
    plotfreq = 1 
    ''')
           
# Main program

#create crystal, relax, run MD/NVE simulation
T_OBJ = 300
initmd(T_OBJ)

makecrystal(4,4,4)

openwindow()

relax_fixbox()
mdpp.cmd('finalcnfile = relaxed.cn writecn')
mdpp.cmd('finalcnfile = relaxed.lammps writeLAMMPS')
mdpp.cmd('finalcnfile = relaxed.cfg writeatomeyecfg')
setup_md()
mdpp.cmd('T_OBJ = {0}'.format(T_OBJ))
mdpp.cmd('initvelocity')
mdpp.cmd('run')

mdpp.cmd('finalcnfile = si100.cn writecn')
mdpp.cmd('finalcnfile = si100.cfg writeatomeyecfg')


# Plot simulation data

from pylab import *

print "Plot simulation data using pylab"

prop = loadtxt('prop.out')
plot( prop[:,0], prop[:,3] )
xlabel("steps")
ylabel("T (K)")
#ion()
show()

# Sleep 
import time
sleep_seconds = 60
print "Python is going to sleep for " + str(sleep_seconds) + " seconds."
time.sleep(sleep_seconds)