薛堪豪的个人博客分享 http://blog.sciencenet.cn/u/bfax 兴趣:凝聚态物理,电化学

博文

VASP结构驰豫过程中检测力和能量收敛的脚本代码(原创)

已有 1120 次阅读 2018-10-9 18:01 |系统分类:科研笔记

用 Fortran 90 编写,代码如下(force.f90):

! This program extracts the force information in each ionic step, during the VASP structural relaxation
! It can be used to monitor force and energy convergence during the relaxation
!
! Compilation:
! > ifort force.f90
! > mv a.out force
!
! Usage:
! In the VASP running path, including OUTCAR
! > force
!
! Written by Kan-Hao Xue, School of Optical and Electronic Information
! Huazhong University of Science and Technology, Wuhan 430074, China
! E-mail: xkh@hust.edu.cn

program force
implicit none

integer :: intAtomCount = 0, intTemp, intStep = 0, istat, i
 real, allocatable :: matForceInfo(:,:)
 real :: maxForce(6)
 character(len=4) index
 character(len=70) line;
 character(len=18) energy;

open(10,file='OUTCAR')

do
 read (10,100, IOSTAT=istat) line
 if (istat < 0) stop 'Cannot find the number of atoms in OUTCAR!'
 if (line(:13) == 'ion position') exit
 end do

do
 read (10,'(A4)', IOSTAT=istat) index
 if (istat < 0) stop 'Cannot find the number of atoms in OUTCAR!'
 if (index ==' LA') exit
 read(index,'(I4)') intTemp
 if (intTemp > intAtomCount) intAtomCount = intTemp
 end do

print '(1X,"Atom count: ", I4)', intAtomCount

allocate(matForceInfo(6, intAtomCount))

print '(/,13X,A,/,13X,A)', 'Max force along each direction', '------------------------------'
 print '(1X,A5,3X,3A11,A13,A18)', 'Index', 'A ', 'B ', 'C ', 'Max force', 'Energy'
 
 do
 read (10,100, IOSTAT=istat) line
 if (istat < 0) exit
 if (line == 'POSITION TOTAL-FORCE (eV/Angst)') then
 read (10,100, IOSTAT=istat) line
 if (istat < 0) stop 'The OUTCAR file is incomplete!'
 read(10,'(6D)') matForceInfo
 matForceInfo = abs(matForceInfo)
 intStep = intStep + 1
 maxForce = maxval(matForceInfo, 2)

do i=1,12
 read (10,100, IOSTAT=istat) line
 if (istat < 0) stop 'The OUTCAR file is incomplete!'
 end do
 
 read (10,200, IOSTAT=istat) energy
 if (istat < 0) stop 'The OUTCAR file is incomplete!'
 
 print '(1X,A4,I3,A1,3F11.5,2X,F11.5,A18)', 'Step', intStep, ':', maxForce(4:6), maxval(maxForce(4:6)), energy
 end if
 end do

close(10)

100 format (1X,A70)
 200 format (64X,A18)
 end


运行实例:

force.fw.png



http://blog.sciencenet.cn/blog-365047-1139878.html

上一篇:Ubuntu 编译和运行程序常见错误处理方法
下一篇:空间群实战分析(二)

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-12-6 13:41

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部