Tool to extract intermediate geometries from OUTCAR

Message

physicalattraction · #1 Post by **physicalattraction** » Mon May 23, 2011 2:42 pm

Is there a tool available which extracts the intermediate geometries from the OUTCAR file and puts them in a multiple structure xyz-file?

#2 Post by **admin** » Mon May 23, 2011 7:57 pm

no, not as far as I know

alex · #3 Post by **alex** » Tue May 24, 2011 11:01 am

#!/usr/bin/awk -f

# (c) by alex, no support
# some issues with many different atom types
# writes an xyz-file on stdout, which is useful for molden or xmol
# it duplicates the cell in a and b direction, not c (useful for surfaces)

BEGIN {

i = 0

# tag for lattice vector reading: read from lattice_read = 1 up to 3
lattice_read = 0
lattice = 0

# tag for atomic position reading
pos_read = 0
position = 0

# tag for opt-cycle
cycle = 0
# tag for finished optcycle (Version 4.5+ only)
optcycle_finished = 0

# tag for ion types
type = 0
number = 0
# initialize
total = 0

}

# READ DATA

{

# collect ions per type from output
# ions per type = 2 4
if ($1 == "ions" && $2 == "per" && $3 == "type") {
for ( j = 1 ; j <= (NF-4) ; j++ ){
ij = 4+j
number++
ionnumber[j]=$ij
# print("number ",number, ionnumber[j])
}
}

# collect ion types from output
if ($1 == "POTCAR:") {
type++
iontype[type]=$3
# print ("type ",type,iontype[type] )
}

# number of dos NEDOS = 301 number of ions NIONS = 68
# read number of ions

if (NF == 12 && $10 == "NIONS") {
nions = $12
}

optcycle_finished++

# FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)
# ---------------------------------------------------
if (NF == 7 && $1 == "FREE" && $2 == "ENERGIE" && $3 == "OF" && $4 == "THE" ) {
optcycle_finished = -2
}

# free energy TOTEN = -593.169883 eV
# read total energy after complete SCF

if (NF == 6 && $3 == "TOTEN" && optcycle_finished == 0) {
cycle++
energy[cycle] = $5
}

# read lattice vector in direct coordinates
# for this we have to read three lines after "direct lattice vectors"

if (lattice_read == 1 ){
a[cycle,lattice,1] = $1
a[cycle,lattice,2] = $2
a[cycle,lattice,3] = $3
# print(a[cycle,lattice,1],a[cycle,lattice,2],a[cycle,lattice,3])

lattice++

if (4==lattice) {
lattice_read = 0
lattice = 0
}
}

if (NF == 6 && $1 == "direct" && $2 == "lattice" && $3 == "vectors"){
lattice_read = 1
lattice = 1
}

# read atomic in cartesian coordinates
# for this we have to read three lines after "direct lattice vectors"

if (pos_read == 1 ){
b[cycle,position-1,1] = $1
b[cycle,position-1,2] = $2
b[cycle,position-1,3] = $3
# print(b[cycle,position,1],b[cycle,position,2],b[cycle,position,3])

position++

if ((position-2) == nions) {
pos_read = 0
position = 0
}
}

if (NF == 3 && $1 == "POSITION" && $2 == "TOTAL-FORCE" && $3 == "(eV/Angst)"){
pos_read = 1
position = 1
}
}

END {

# generate type/number pairs (easier to handle during printing)
for ( k = 1 ; k <= number ; k++ ){
for ( j = 1 ; j <= ionnumber[k] ; j++ ){
total++
ion[total] = iontype[k]
# print(total, ion[total], k,j,iontype[k])
}
}

# print all information per cycle

for (cyc = 1 ; cyc <= cycle-1 ; cyc++) {

# print number of ions
print nions

# print cycle & energy
printf("ene:%8i%15.6f%s\n", cyc, energy[cyc], " eV")

# print ionposition in cartesian coordinates
for ( j = 1 ; j <= nions ; j++ ){
printf("%s %10.5f%10.5f%10.5f\n", ion[j], b[cyc,j,1], b[cyc,j,2], b[cyc,j,3])
}
}

}

physicalattraction · #4 Post by **physicalattraction** » Wed May 25, 2011 7:19 pm

Thank you for your bash script. I have now written the following fortran program, which extracts the positions I need.

Code: Select all

c??Reads?VASP?OUTCAR?file
??????integer,?parameter::?natom?=?27
??????character(len=70)?::?label,?num,?fname
??????character(len=2)??::?symb(natom)

??????real*8????????????::?vec(3),?cell(3,3),?sh(3)

c-----------------------------------------------------------------------

??????symb(:)(:)?=?'?'
??????do?ia?=?1,?6
???????symb(ia)=?'C'
??????enddo
??????do?ia?=?7,?11
???????symb(ia)=?'H'
??????enddo
??????symb(12)=?'S'
??????do?ia?=?13,?natom
???????symb(ia)=?'Ag'
??????enddo

??????nargs?=?iargc()
??????if?(?nargs.lt.1?)?stop?'Usage:?file_names?XV2MD?shfit(1:3)?fac'
??????call?getarg(1,fname)

??????ncel1=?0;?ncel2=0;?ncel3=0;?fac=1.0
??????if?(?nargs.gt.1?)?then
???????call?getarg(2,num);?read(num,*)?ncel1
??????endif
??????if?(?nargs.gt.2?)?then
???????call?getarg(3,num);?read(num,*)?ncel2
??????endif
??????if?(?nargs.gt.3?)?then
???????call?getarg(4,num);?read(num,*)?ncel3
??????endif
??????write(6,*)?'ncel=',ncel1,ncel2,ncel3
??????if?(?nargs.gt.4?)?then
???????call?getarg(5,num);?read(num,*)?fac
??????endif
??????write(6,*)?'fac=',fac
??????open(1,file='OUTCAR',status='old',err=3)
??????open(2,file=trim(fname)//'_MD.xyz',status='unknown')
??????lcell?=?.false.
??????do
???????read(1,'(a)',end=10)?label

c??????read?cell?
???????if?(?.not.?lcell?)?then
????????if?(?index(label,'direct?lattice?vectors').ne.0?)?then
?????????do?i?=?1?,?3
??????????read(1,*)?cell(:,i)
?????????enddo
?????????lcell?=?.true.
????????endif
???????endif

c??????read?positions?
???????if?(?index(label,'POSITION').ne.0?)?then
????????read(1,*)
????????natoms?=?(2*ncel1+1)*(2*ncel2+1)*(2*ncel3+1)*natom
????????write(2,'(i6/)')?natoms
????????do?ia?=?1,?natom
?????????read(1,*)?vec
?????????if?(?ia.gt.12?)?vec?=?vec?-?cell(:,3)
?????????do?icel3=?-ncel3,?ncel3
??????????do?icel2=?-ncel2,?ncel2
???????????do?icel1=?-ncel1,?ncel1
????????????sh?=?icel1*cell(:,1)??+?icel2*cell(:,2)?+?icel3*cell(:,3)
????????????write(2,'(1x,a,3f10.5)')?symb(ia),vec(:)+sh
???????????enddo
??????????enddo
?????????enddo
????????enddo
???????endif
??????enddo
10????close(1)
??????close(2)

??????stop
3?????stop?'could?not?open?OUTCAR?file'
??????end

physicalattraction · #5 Post by **physicalattraction** » Wed May 25, 2011 7:21 pm

Hmmm, in the preview, it didn't show every space character as a question mark. There is probably a bug in the forum software. I can't edit the message above anymore, but if someone uses it: every question mark is a space.