# -------------------------------------------------------------
# Custom Analysis
# -------------------------------------------------------------
md_trajectory = nlread('polycrystal_creep.nc')[-1]
# Extract the times, cell vectors, potential energies, and pressure tensors
times        = md_trajectory.times().inUnitsOf(picoSecond)
cell_vectors = md_trajectory.cells()
e_pot        = md_trajectory.potentialEnergies().inUnitsOf(eV)
pressure       = md_trajectory.pressureTensors().inUnitsOf(GPa)
# Convert the tensor elements into lists
strain_xx = cell_vectors[:, 0, 0]/cell_vectors[0, 0, 0] - 1.0
strain_zz = cell_vectors[:, 2, 2]/cell_vectors[0, 2, 2] - 1.0
pressure_zz = pressure[:, 2, 2]
# Plot the data usig pylab
import pylab
# Plot the strain
pylab.figure()
pylab.plot(times,strain_zz,'b')
pylab.plot(times,strain_xx,'g')
pylab.title('Creep strain in Copper Polycrystal at 500K and a load of 1.5 GPa')
pylab.xlabel('MD time (ps)')
pylab.ylabel('Strain')
pylab.grid(True)
pylab.savefig("CreepStrain.png")
# Plot the pressure
pylab.figure()
pylab.plot(times,pressure_zz,'b')
pylab.title('Stress in Copper Polycrystal')
pylab.xlabel('MD time (ps)')
pylab.ylabel('Stress (GPa,bar)')
pylab.grid(True)
pylab.savefig("CreepStress.png")
# Plot the potential energy
pylab.figure()
pylab.plot(times,e_pot)
pylab.title('Potential energy during creep in copper polycrystal')
pylab.xlabel('MD time (ps)')
pylab.ylabel('Potential energy (eV)')
pylab.grid(True)
pylab.savefig("CreepEPot.png" )