Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Save stress tensor in extXYZ file #379

Merged
merged 2 commits into from
Nov 5, 2024
Merged

Save stress tensor in extXYZ file #379

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f873ac9
Add stress to extXYZ file.
AugustinLu Nov 4, 2024
46899ed
Modify unit of output stress from eV to eV/ang**3 by division by cell…
AugustinLu Nov 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 14 additions & 12 deletions src/control.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,8 @@ module control
!! coordinates.
!! 2022/09/16 16:51 dave
!! Added SQNM options for ion/cell optimisation (only partly complete)
!! 2024/11/04 Augustin Lu
!! Added stress to write_extxyz call
!! SOURCE
!!
subroutine control_run(fixed_potential, vary_mu, total_energy)
Expand All @@ -147,7 +149,7 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
use dimens, only: r_core_squared, r_h
use GenComms, only: my_barrier, cq_abort, cq_warn
use pseudopotential_data, only: set_pseudopotential
use force_module, only: tot_force
use force_module, only: tot_force, stress
use minimise, only: get_E_and_F
use global_module, only: runtype, flag_self_consistent, &
flag_out_wf, flag_write_DOS, wf_self_con, &
Expand Down Expand Up @@ -188,7 +190,7 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
endif
call get_E_and_F(fixed_potential, vary_mu, total_energy,&
flag_ff, flag_wf, level=backtrace_level)
if (flag_write_extxyz) call write_extxyz('trajectory.xyz', total_energy, tot_force)
if (flag_write_extxyz) call write_extxyz('trajectory.xyz', total_energy, tot_force,stress)
!
else if ( leqi(runtype, 'cg') ) then
if (flag_opt_cell) then
Expand Down Expand Up @@ -320,7 +322,7 @@ subroutine cg_run(fixed_potential, vary_mu, total_energy)
safemin2, cg_line_min, safe, adapt_backtrack, backtrack
use GenComms, only: gsum, myid, inode, ionode
use GenBlas, only: dot
use force_module, only: tot_force
use force_module, only: tot_force, stress
use io_module, only: write_atomic_positions, pdb_template, &
check_stop, write_xsf, write_extxyz, &
print_atomic_positions, return_prefix
Expand Down Expand Up @@ -496,7 +498,7 @@ subroutine cg_run(fixed_potential, vary_mu, total_energy)
end if
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Analyse forces
g0 = dot(length, tot_force, 1, tot_force, 1)
call get_maxf(max)
Expand Down Expand Up @@ -1558,7 +1560,7 @@ subroutine lbfgs(fixed_potential, vary_mu, total_energy)
use move_atoms, only: updateL, updateLorK, updateSFcoeff
use GenComms, only: gsum, myid, inode, ionode, gcopy, my_barrier
use GenBlas, only: dot
use force_module, only: tot_force
use force_module, only: tot_force, stress
use io_module, only: write_atomic_positions, pdb_template, &
check_stop, write_xsf, write_extxyz
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl
Expand Down Expand Up @@ -1718,7 +1720,7 @@ subroutine lbfgs(fixed_potential, vary_mu, total_energy)
! Add call to write_atomic_positions and write_xsf (2020/01/17: smujahed)
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Build q
do i=iter, iter_low, -1
! Indexing
Expand Down Expand Up @@ -1848,7 +1850,7 @@ subroutine sqnm(fixed_potential, vary_mu, total_energy)
use move_atoms, only: single_step, backtrack_linemin
use GenComms, only: myid, inode, ionode
use GenBlas, only: dot, syev
use force_module, only: tot_force
use force_module, only: tot_force, stress
use io_module, only: write_atomic_positions, pdb_template, write_extxyz, &
check_stop, write_xsf, print_atomic_positions, return_prefix
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl
Expand Down Expand Up @@ -2028,7 +2030,7 @@ subroutine sqnm(fixed_potential, vary_mu, total_energy)
! Add call to write_atomic_positions and write_xsf (2020/01/17: smujahed)
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Build significant subspace
Sij = zero
omega = zero
Expand Down Expand Up @@ -2438,7 +2440,7 @@ subroutine cell_sqnm(fixed_potential, vary_mu, total_energy)
! Add call to write_atomic_positions and write_xsf (2020/01/17: smujahed)
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Build significant subspace
Sij = zero
omega = zero
Expand Down Expand Up @@ -2888,7 +2890,7 @@ subroutine full_sqnm(fixed_potential, vary_mu, total_energy)
! Add call to write_atomic_positions and write_xsf (2020/01/17: smujahed)
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Build significant subspace
Sij = zero
omega = zero
Expand Down Expand Up @@ -3288,7 +3290,7 @@ subroutine cell_cg_run(fixed_potential, vary_mu, total_energy)
rcellx, d_units(dist_units), rcelly, d_units(dist_units), rcellz, d_units(dist_units)
end if
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)
! Analyse Stresses and energies
dH = enthalpy1 - enthalpy0
volume = rcellx*rcelly*rcellz
Expand Down Expand Up @@ -4356,7 +4358,7 @@ subroutine full_cg_run_single_vector(fixed_potential, vary_mu, total_energy)
end if
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force)
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) call write_extxyz('trajectory.xyz', energy1, tot_force, stress)

! Analyse forces and stress
g0 = dot(length-3, tot_force, 1, tot_force, 1)
Expand Down
25 changes: 22 additions & 3 deletions src/io_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2870,10 +2870,11 @@ end subroutine write_xsf
!! CREATION DATE
!! 2021/10/18
!! MODIFICATION HISTORY
!!
!! 2024/11/04 Augustin Lu
!! Add stress tensor
!! SOURCE
!!
subroutine write_extxyz(filename, energy0, atom_force)
subroutine write_extxyz(filename, energy0, atom_force, stress_tensor)

use datatypes
use timer_module
Expand All @@ -2889,14 +2890,17 @@ subroutine write_extxyz(filename, energy0, atom_force)
character(len=*) :: filename
real(double) :: energy0
real(double), dimension(3,ni_in_cell) :: atom_force
real(double), dimension(3,3) :: stress_tensor

! Local variables
integer :: lun, i, j, title_length
character(len=2) :: atom_name
character(len=432) :: comment
character(len=45) :: vec_a, vec_b, vec_c, energy_str
character(len=80) :: titles_xyz
real(double) :: for_conv_loc, en_conv_loc, dist_conv_loc
character(len=135) :: stress_str

real(double) :: for_conv_loc, en_conv_loc, dist_conv_loc, volume

if(inode==ionode) then
if (iprint_init>2) write(io_lun, &
Expand Down Expand Up @@ -2930,6 +2934,21 @@ subroutine write_extxyz(filename, energy0, atom_force)
comment=TRIM(comment)//' Properties=species:S:1:pos:R:3:forces:R:3 potential_energy='
write(energy_str,'(f0.8)') energy0 * en_conv_loc
comment = TRIM(comment)//TRIM(energy_str)//' pbc="T T T" '

volume = r_super_x*r_super_y*r_super_z*BohrToAng**3

! Convert each row to string
write(stress_str(1:45), '(3f15.8)') stress_tensor(1,1)*HaToeV/volume,&
stress_tensor(1,2)*HaToeV/volume,&
stress_tensor(1,3)*HaToeV/volume
write(stress_str(46:90), '(3f15.8)') stress_tensor(2,1)*HaToeV/volume,&
stress_tensor(2,2)*HaToeV/volume,&
stress_tensor(2,3)*HaToeV/volume
write(stress_str(91:135), '(3f15.8)') stress_tensor(3,1)*HaToeV/volume,&
stress_tensor(3,2)*HaToeV/volume,&
stress_tensor(3,3)*HaToeV/volume
comment = TRIM(comment)//' stress=" '//TRIM(stress_str)//'" '

write(lun,'(a)') TRIM(comment)

do i=1,ni_in_cell
Expand Down
10 changes: 6 additions & 4 deletions src/md_misc_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@ end subroutine init_md
!! 2020/10/07 tsuyoshi
!! added deallocation of atom_vels
!! SOURCE
!!
!!
subroutine end_md(th, baro)

use GenComms, only: inode, ionode
Expand Down Expand Up @@ -313,7 +313,7 @@ end subroutine end_md
!! CREATION DATE
!! 2018/04/23 11:49
!! SOURCE
!!
!!
subroutine integrate_pt(baro, thermo, mdl, velocity, second_call)

use numbers
Expand Down Expand Up @@ -551,8 +551,10 @@ end subroutine get_heat_flux
!! MODIFIED
!! 2021/10/19 Jianbo Lin
!! Added call for extended XYZ output (includes forces)
!! 2024/11/04 Augustin Lu
!! Added stress in call for extended XYZ output
!! SOURCE
!!
!!
subroutine write_md_data(iter, thermo, baro, mdl, nequil, MDfreq, XSFfreq, XYZfreq)

use GenComms, only: inode, ionode
Expand Down Expand Up @@ -582,7 +584,7 @@ subroutine write_md_data(iter, thermo, baro, mdl, nequil, MDfreq, XSFfreq, XYZfr
call write_xsf(md_trajectory_file, iter)
end if
if (flag_write_extxyz .and. mod(iter,XYZfreq) == 0) then
call write_extxyz('trajectory.xyz', mdl%dft_total_energy, mdl%atom_force)
call write_extxyz('trajectory.xyz', mdl%dft_total_energy, mdl%atom_force, mdl%stress)
end if
if (flag_heat_flux) call mdl%dump_heat_flux(md_heat_flux_file)
if (mod(iter, MDfreq) == 0) then
Expand Down