-
Notifications
You must be signed in to change notification settings - Fork 0
/
chgcar.f
162 lines (141 loc) · 3.4 KB
/
chgcar.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
program chgcar
use ifport
use prec
use mio
use mchgcar
implicit none
! input and output units
integer,dimension(3) :: iou = (/50,51,52/)
character(len=32),dimension(3) :: ionam
character(len=32767),dimension(2) :: line
character(len=32767) :: buf
real(q),dimension(2) :: c
real(q),dimension(2,3) :: pos
real(q),dimension(3,5) :: v
integer,dimension(2,4) :: ng
integer :: i,j,ni,lm,lr,lc,ioerr
! timing
real :: totalt
real,dimension(2) :: telapsed
! structures
type(t_poscar),dimension(2) :: car
print*,'Warning: Spin density is not implemented!'
lparam: do i=1,3
if(i.lt.3) read(*,fmt="(f3.6)"),c(i)
read(*,fmt="(a)"),ionam(i)
end do lparam
!
lopen: do i=1,3
if(i.eq.3) then
call wopen(iou(i),ionam(i),ioerr)
else
call ropen(iou(i),ionam(i),ioerr)
end if
if(ioerr/=0) then
print*,"Open error: ",ionam(i),ioerr
stop
end if
end do lopen
! read header
lrdhdr: do i=1,2
call chgrdhdr(iou(i),car(i),ioerr)
if(ioerr/=0) then
print*,"Header error: ",ionam(i),ioerr
stop
end if
end do lrdhdr
! check 1
if(car(1)%nion.ne.car(2)%nion) then
print*,"Warning: Atom numbers are not equal!",car(1)%nion,car(2)%nion
! stop
end if
print*
! write output header
call chgwrhdr(iou(3),car(1),ioerr)
if(ioerr/=0) then
print*,"Header error: ",ionam(3),ioerr
stop
end if
! dump coords
lcoord1: do ni=1,car(1)%nion
read(iou(1),*,iostat=ioerr) (pos(1,i),i=1,3)
if(ioerr/=0) stop
write(unit=iou(3),fmt="(3f12.6)",iostat=ioerr) (pos(1,i),i=1,3)
if(ioerr/=0) stop
end do lcoord1
! pad read 2 input
lcoord2: do ni=1,car(2)%nion
read(iou(2),*,iostat=ioerr) (pos(2,i),i=1,3)
if(ioerr/=0) stop
end do lcoord2
! dump padding
write(unit=iou(3),fmt="(a)",iostat=ioerr) ''
if(ioerr/=0) stop
! grid dim
lng: do i=1,2
read(iou(i),*,iostat=ioerr) (ng(i,j),j=1,3)
if(ioerr/=0) stop
print "(f6.3,x,20a)",c(i),ionam(i)
print "(3i5)",(ng(i,j),j=1,3)
print*
end do lng
ng(1,4) = 1
ng(2,4) = 1
! check
do j=1,3
if(ng(1,j).ne.ng(2,j)) then
print*,'Grid dimension error'
stop
end if
do i=1,2
ng(i,4) = ng(i,4) * ng(i,j)
end do
end do
write(unit=iou(3),fmt="(3i8)",iostat=ioerr) (ng(1,i),i=1,3)
if(ioerr/=0) stop
! grid data
lm=floor(ng(1,4)/5.0)
lr=mod(ng(1,4),5)
! print*,ng(1,4),lm*5,lm,lr
print*,"Processing..."
lgrid: do lc=1,lm
call vline(iou,ioerr,c,v,5)
end do lgrid
if(lr.gt.0) then
call vline(iou,ioerr,c,v,lr)
end if
!
lclose: do i=1,3
close(unit=iou(i))
end do lclose
totalt=etime(telapsed)
print "(f12.3,x,6a)",totalt,"sec"
contains
subroutine vcomp(c,v,im)
real(q),dimension(2) :: c
real(q),dimension(3,5) :: v
integer :: i,im
do i=1,im
v(3,i)=c(1)*v(1,i)+c(2)*v(2,i)
end do
end subroutine
subroutine vline(iou,ioerr,c,v,jm)
real(q),dimension(2) :: c
real(q),dimension(3,5) :: v
integer :: i,j,jm,ioerr
integer,dimension(3) :: iou
do i=1,2
read(iou(i),*,iostat=ioerr) (v(i,j),j=1,jm)
if(ioerr/=0) then
print*,"Input error: ",i,ioerr
stop
end if
end do
call vcomp(c,v,jm)
write(iou(3),fmt="(5E18.9)",iostat=ioerr) (v(3,j),j=1,jm)
if(ioerr/=0) then
print*,"Output error: ",i,ioerr
stop
end if
end subroutine
end program