-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathunwrap_PBC_cell.f90
131 lines (115 loc) · 4.03 KB
/
unwrap_PBC_cell.f90
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
! program unwrap_PBC.f
!******************************************************************
! Unwrap periodically wrapped coordinates from an xyz file
! Output the unwrapped coords to a new xyz file
!
! NOTE: This algorithm assumes a reasonably small stepsize
! with respect to the lattice constants. Is partially
! taken from F.26 from Allen-Tildesley.
! I. Tamblyn modified B.Boates code to allow for variable cell
! size over the course of a simulation (for use with NPT / shock
! simulations
!******************************************************************
implicit none
! Declare all required variables
! n = max # of timesteps, m = max # of atoms
integer n, m, i, j, k
parameter (n=100000,m=512)
integer natoms(n), timesteps(n)
character(2) typat(n,m)
character(128) fin
character(8) dummy1, dummy2, dummy3
real(4) Rx(n,m),Ry(n,m),Rz(n,m)
real(4) dx(n,m),dy(n,m),dz(n,m),dRx,dRy,dRz
real(4) a, ax, ay, az
! Ask user for required input, which can be given via
! an input file (still echoes the requests though...)
write(*,*) "Name of xyz file with wrapped coordinates"
read(*,*) fin
! write(*,*) "Lattice constants (ax,ay,az)"
! read(*,*) ax, ay, az
! Assume cubic cell:
! ax = a
! ay = a
! az = a
! Open the chosen input file in "unit" 1
! Which must be simply 3 columns of coordinates
open(1,file=fin,status='old',ERR=90)
! Open an output file of the chosen name in "unit" 2
open(2,file='unwrapped.xyz')
! Read in the wrapped xyz file data
! Convert to reduced coordinates and shift origin
do i=1,n
read(1,*,END=100) natoms(i)
read(1,*) ax, ay, az, dummy1, dummy2, dummy3, timesteps(i)
do j=1,natoms(i)
read(1,*) typat(i,j),Rx(i,j),Ry(i,j),Rz(i,j)
Rx(i,j) = (Rx(i,j) / ax) + 0.5
Ry(i,j) = (Ry(i,j) / ay) + 0.5
Rz(i,j) = (Rz(i,j) / az) + 0.5
enddo
enddo
100 continue
! Find all dx,dy,dz values from the wrapped data
do j=1,natoms(1)
dx(1,j) = 0.0
dy(1,j) = 0.0
dz(1,j) = 0.0
enddo
do i=2,n
do j=1,natoms(i)
dRx = Rx(i,j) - Rx(i-1,j)
dRy = Ry(i,j) - Ry(i-1,j)
dRz = Rz(i,j) - Rz(i-1,j)
if (dRx.ge.ax/2) then
print *, "WARNING: dRx has exceeded a/2"
endif
if (dRy.ge.ay/2) then
print *, "WARNING: dRy has exceeded a/2"
endif
if (dRz.ge.az/2) then
print *, "WARNING: dRz has exceeded a/2"
endif
dx(i,j) = dx(i-1,j) + (dRx - anint(dRx))
dy(i,j) = dy(i-1,j) + (dRy - anint(dRy))
dz(i,j) = dz(i-1,j) + (dRz - anint(dRz))
enddo
enddo
!$omp parallel do
!$omp& default(private)
!$omp& shared(natoms,timesteps)
! Unwrap the data
do i=1,n
if (natoms(i).eq.0) then
goto 80
endif
do j=1,natoms(i)
! Calculate the unwrapped coordinate
Rx(i,j) = Rx(1,j) + dx(i,j)
Ry(i,j) = Ry(1,j) + dy(i,j)
Rz(i,j) = Rz(1,j) + dz(i,j)
enddo
enddo
80 continue
!$omp end parallel do
! Write to file
do i=1,n
if (natoms(i).eq.0) then
goto 70
endif
write(2,*) natoms(i)
write(2,*) ax, ay, az
do j=1,natoms(i)
! Convert back to Cartesian coordinates (re-shift origin)
Rx(i,j) = (Rx(i,j) - 0.5) * ax
Ry(i,j) = (Ry(i,j) - 0.5) * ay
Rz(i,j) = (Rz(i,j) - 0.5) * az
! Write unwrapped coordinates to file
write(2,*) typat(i,j), Rx(i,j), Ry(i,j), Rz(i,j)
enddo
enddo
70 continue
close(1)
close(2)
90 continue
END