Examples and docs

timestep	0.0001

# compute and output options

compute out_mag    all compute/spin

compute out_pe     all pe

#!/usr/bin/env python

# This file reads in the file log.lammps generated by the script ELASTIC/in.elastic

# It prints out the 6x6 tensor of elastic constants Cij 

# followed by the 6x6 tensor of compliance constants Sij

# It uses the same conventions as described in:

# Sprik, Impey and Klein PRB (1984).  

# The units of Cij are whatever was used in log.lammps (usually GPa)

# The units of Sij are the inverse of that (usually 1/GPa)

from numpy import zeros

from numpy.linalg import inv

# define logfile layout

nvals = 21

valpos = 4

valstr = '\nElastic Constant C'

# define order of Cij in logfile

cindices = [0]*nvals

cindices[0] = (0,0)

cindices[1] = (1,1)

cindices[2] = (2,2)

cindices[3] = (0,1)

cindices[4] = (0,2)

cindices[5] = (1,2)

cindices[6] = (3,3)

cindices[7] = (4,4)

cindices[8] = (5,5)

cindices[9] = (0,3)

cindices[10] = (0,4)

cindices[11] = (0,5)

cindices[12] = (1,3)

cindices[13] = (1,4)

cindices[14] = (1,5)

cindices[15] = (2,3)

cindices[16] = (2,4)

cindices[17] = (2,5)

cindices[18] = (3,4)

cindices[19] = (3,5)

cindices[20] = (4,5)

# open logfile

logfile = open("log.lammps",'r')

txt = logfile.read()

# search for 21 elastic constants

c = zeros((6,6))

s2 = 0

for ival in range(nvals):

    s1 = txt.find(valstr,s2)

    if (s1 == -1):

        print "Failed to find elastic constants in log file"

        exit(1)

    s1 += 1

    s2 = txt.find("\n",s1)

    line = txt[s1:s2]

#    print line

    words = line.split()

    (i1,i2) = cindices[ival]

    c[i1,i2] = float(words[valpos])

    c[i2,i1] = c[i1,i2]

print "C tensor [GPa]"

for i in range(6):

    for j in range(6):

        print "%10.8g " % c[i][j],

    print

# apply factor of 2 to columns of off-diagonal elements

for i in range(6):

    for j in range(3,6):

        c[i][j] *= 2.0

s = inv(c)

# apply factor of 1/2 to columns of off-diagonal elements

for i in range(6):

    for j in range(3,6):

        s[i][j] *= 0.5

print "S tensor [1/GPa]"

for i in range(6):

    for j in range(6):

        print "%10.8g " % s[i][j],

    print

# NOTE: This script should not need to be

# modified. See in.elastic for more info.

#

# Find which reference length to use

if "${dir} == 1" then &

   "variable len0 equal ${lx0}" 

if "${dir} == 2" then &

   "variable len0 equal ${ly0}" 

if "${dir} == 3" then &

   "variable len0 equal ${lz0}" 

if "${dir} == 4" then &

   "variable len0 equal ${lz0}" 

if "${dir} == 5" then &

   "variable len0 equal ${lz0}" 

if "${dir} == 6" then &

   "variable len0 equal ${ly0}" 

# Reset box and simulation parameters

clear

box tilt large

read_restart restart.equil

include potential.mod

# Negative deformation

variable delta equal -${up}*${len0}

variable deltaxy equal -${up}*xy

variable deltaxz equal -${up}*xz

variable deltayz equal -${up}*yz

if "${dir} == 1" then &

   "change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box"

if "${dir} == 2" then &

   "change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box"

if "${dir} == 3" then &

   "change_box all z delta 0 ${delta} remap units box"

if "${dir} == 4" then &

   "change_box all yz delta ${delta} remap units box"

if "${dir} == 5" then &

   "change_box all xz delta ${delta} remap units box"

if "${dir} == 6" then &

   "change_box all xy delta ${delta} remap units box"

# Relax atoms positions

minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Obtain new stress tensor

variable tmp equal pxx

variable pxx1 equal ${tmp}

variable tmp equal pyy

variable pyy1 equal ${tmp}

variable tmp equal pzz

variable pzz1 equal ${tmp}

variable tmp equal pxy

variable pxy1 equal ${tmp}

variable tmp equal pxz

variable pxz1 equal ${tmp}

variable tmp equal pyz

variable pyz1 equal ${tmp}

# Compute elastic constant from pressure tensor

variable C1neg equal ${d1}

variable C2neg equal ${d2}

variable C3neg equal ${d3}

variable C4neg equal ${d4}

variable C5neg equal ${d5}

variable C6neg equal ${d6}

# Reset box and simulation parameters

clear

box tilt large

read_restart restart.equil

include potential.mod

# Positive deformation

variable delta equal ${up}*${len0}

variable deltaxy equal ${up}*xy

variable deltaxz equal ${up}*xz

variable deltayz equal ${up}*yz

if "${dir} == 1" then &

   "change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box"

if "${dir} == 2" then &

   "change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box"

if "${dir} == 3" then &

   "change_box all z delta 0 ${delta} remap units box"

if "${dir} == 4" then &

   "change_box all yz delta ${delta} remap units box"

if "${dir} == 5" then &

   "change_box all xz delta ${delta} remap units box"

if "${dir} == 6" then &

   "change_box all xy delta ${delta} remap units box"

# Relax atoms positions

minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Obtain new stress tensor

variable tmp equal pe

variable e1 equal ${tmp}

variable tmp equal press

variable p1 equal ${tmp}

variable tmp equal pxx

variable pxx1 equal ${tmp}

variable tmp equal pyy

variable pyy1 equal ${tmp}

variable tmp equal pzz

variable pzz1 equal ${tmp}

variable tmp equal pxy

variable pxy1 equal ${tmp}

variable tmp equal pxz

variable pxz1 equal ${tmp}

variable tmp equal pyz

variable pyz1 equal ${tmp}

# Compute elastic constant from pressure tensor

variable C1pos equal ${d1}

variable C2pos equal ${d2}

variable C3pos equal ${d3}

variable C4pos equal ${d4}

variable C5pos equal ${d5}

variable C6pos equal ${d6}

# Combine positive and negative 

variable C1${dir} equal 0.5*(${C1neg}+${C1pos})

variable C2${dir} equal 0.5*(${C2neg}+${C2pos})

variable C3${dir} equal 0.5*(${C3neg}+${C3pos})

variable C4${dir} equal 0.5*(${C4neg}+${C4pos})

variable C5${dir} equal 0.5*(${C5neg}+${C5pos})

variable C6${dir} equal 0.5*(${C6neg}+${C6pos})

# Delete dir to make sure it is not reused

variable dir delete

# fcc cobalt in a 3d periodic box

clear 

units metal

atom_style spin

dimension 3

boundary p p p

# check why?

atom_modify map array 

lattice fcc 3.54

region box block 0.0 5.0 0.0 5.0 0.0 5.0

create_box 1 box

create_atoms 1 box

# setting mass, mag. moments, and interactions for cobalt

mass		1 58.93

set group all spin/random 31 1.72

velocity all create 100 4928459 rot yes dist gaussian

#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0

pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0

pair_coeff * * eam/alloy ../examples/SPIN/cobalt/Co_PurjaPun_2012.eam.alloy Co

pair_coeff * * pair/spin/exchange exchange 4.0 0.0446928 0.003496 1.4885 

#pair_coeff * * pair/spin/soc/neel neel 4.0 0.003330282 0.864159 2.13731 

neighbor 0.1 bin

neigh_modify every 10 check yes delay 20

fix 1 all force/spin zeeman 0.1 0.0 0.0 1.0

fix 2 all langevin/spin 0.0 0.1 21

fix 3 all integration/spin lattice yes

timestep	0.0001

compute out_mag    all compute/spin

compute out_pe     all pe

compute out_ke     all ke

compute out_temp   all temp

variable magz      equal c_out_mag[3]

variable magnorm   equal c_out_mag[4]

variable emag      equal c_out_mag[5]

variable tmag      equal c_out_mag[6]

variable mag_force equal f_1

thermo_style    custom step time v_magnorm v_emag temp etotal

thermo          10

#dump 1 all custom 50 dump_cobalt.lammpstrj type x y z spx spy spz

run 1000