Statistics
| Revision:

root / src / CheckPeriodicBound.f90 @ 10

History | View | Annotate | Download (6.3 kB)

1
SUBROUTINE CheckPeriodicBound
2
  ! This subroutine is only for periodic systems
3
  ! it checks that the atoms do not change their
4
  ! positions from one side to the other during the simulations.
5
  ! This would ruin the interpolation.
6

    
7

    
8
  Use Path_module
9
  Use Io_module
10

    
11
  IMPLICIT NONE
12

    
13

    
14
  INTEGER(KINT) :: IGeom, I,J,K
15
  REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat)
16

    
17
  REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt
18
  REAL(KREAL) :: NormA, NormB,Xrel,Yrel,ZRel,NormC
19
  REAL(KREAL) :: Prod,LatrA(3),LatrB(3)
20
  REAL(KREAL) :: V(3)
21

    
22
  LOGICAL :: Debug
23

    
24
  INTERFACE
25
     function valid(string) result (isValid)
26
       CHARACTER(*), intent(in) :: string
27
       logical                  :: isValid
28
     END function VALID
29

    
30

    
31
    SUBROUTINE die(routine, msg, file, line, unit)
32

    
33
      Use VarTypes
34
      Use io_module
35

    
36
      implicit none
37
      character(len=*), intent(in)           :: routine, msg
38
      character(len=*), intent(in), optional :: file
39
      integer(KINT), intent(in), optional      :: line, unit
40

    
41
    END SUBROUTINE die
42

    
43
  END INTERFACE
44

    
45
  debug=valid('checkperiodicbound')
46
  if (debug) call header("Entering CheckPeriodicBound")
47

    
48
  ALLOCATE(GeomRef(3,Nat), GeomTmp(3,Nat))
49

    
50
  ! The test is easy in DIRECT space, so we will convert the geometries
51
  ! into this space. 
52

    
53
  SELECT CASE (Iper)
54
  CASE(1)
55
     NormA=DOT_PRODUCT(Lat_a,Lat_a)
56
     DO IGeom=2,NGeomI
57
        DO I=1,NAt
58
           DO J=1,3
59
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
60
           END DO
61
           XRel=DOT_PRODUCT(Lat_a,v)/NormA
62
           if (abs(Xrel).GE.BoxTol) THEN
63
              If (debug) THEN
64
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
65
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
66
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
67
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)
68
              END IF
69
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)
70
           END IF
71
        END DO
72
     END DO
73
  CASE (2)
74
     NormA=DOT_PRODUCT(Lat_a,Lat_a)
75
     NormB=DOT_PRODUCT(Lat_b,Lat_b)
76
     Prod=DOT_PRODUCT(Lat_a,Lat_b)
77
     LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)
78
     LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)
79

    
80
     DO IGeom=2,NGeomI
81
        DO I=1,NAt
82
           DO J=1,3
83
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
84
           END DO
85
           XRel=DOT_PRODUCT(Latra,v)
86
           YRel=DOT_PRODUCT(Latrb,v)
87
           if (abs(Xrel).GE.BoxTol) THEN
88
              If (debug) THEN
89
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
90
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
91
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
92
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)
93
              END IF
94
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)
95
           END IF
96
           if (abs(Yrel).GE.BoxTol) THEN
97
              If (debug) THEN
98
                 WRITE(*,'("Atom ",I5," moved because YRel=",F10.6)') I,YRel
99
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
100
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
101
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)
102
              END IF
103
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)
104
           END IF
105
        END DO
106
     END DO
107

    
108
  CASE(3)
109

    
110

    
111

    
112
!      if (debug) THEN
113
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_a:",Lat_a
114
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_b:",Lat_b
115
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_c:",Lat_c
116
!         WRITE(*,'(1X,A,3(1X,F10.6))') "BocTol:",BoxTol
117

    
118
!      NormA=DOT_PRODUCT(Lat_a,Lat_a)
119
!      NormB=DOT_PRODUCT(Lat_b,Lat_b)
120
!      Prod=DOT_PRODUCT(Lat_a,Lat_b)
121
!      LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)
122
!      LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)
123
! !     NormA=DOT_PRODUCT(Latra,Latra)
124
! !     NormB=DOT_PRODUCT(Latrb,Latrb)
125

    
126
!         DO IGeom=2,NGeomI
127
!            DO I=1,NAt
128
!               DO J=1,3
129
!                  v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
130
!               END DO
131
!               XRel=DOT_PRODUCT(Latra,v)
132
!               YRel=DOT_PRODUCT(Latrb,v)
133
!               If (abs(XRel).GT.BoxTol)  WRITE(*,*) "Xrel too big:",IGeom,I,Xrel
134
!               If (abs(YRel).GT.BoxTol)  WRITE(*,*) "Yrel too big:",IGeom,I,Yrel
135
! !               WRITE(*,*) "Xrel :",IGeom,I,Xrel
136
! !               WRITE(*,*) "Yrel :",IGeom,I,Yrel
137
! !               WRITE(*,*) "Zrel :",IGeom,I,Zrel
138

    
139
!            END DO
140
!         END DO
141
!      END IF
142

    
143
     Latrloc(1:3,1)=Lat_a
144
     Latrloc(1:3,2)=Lat_b
145
     Latrloc(1:3,3)=Lat_c
146
     Bloc=1.
147
     ! TODO: Replace by LAPACK subroutine
148
     CALL Gaussj(Latrloc,3,3,Bloc,1,1)
149

    
150
     ! We put first geometry as the reference
151
     GeomRef=0.
152
     DO I=1,Nat
153
        DO K=1,3
154
           DO J=1,3
155
              GeomRef(K,I)=GeomRef(K,I)+XyzGeomI(1,J,I)*LatrLoc(K,J)
156
           END DO
157
        END DO
158
     END DO
159

    
160
     DO IGeom=2,NGeomI
161
        GeomTmp=0.
162
        DO I=1,Nat
163
           DO K=1,3
164
              DO J=1,3
165
                 GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J)
166
              END DO
167
              if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN
168
                 If (debug) & 
169
                      WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &
170
                      K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),& 
171
                      GeomRef(K,I),GeomTmp(K,I), &
172
                      GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))
173

    
174
                 GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))
175
              END IF
176
           END DO
177
           Xt=GeomTmp(1,I)
178
           Yt=GeomTmp(2,I)
179
           Zt=GeomTmp(3,I)
180
           XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)
181
           XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)
182
           XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)
183
        END DO
184
     END DO
185

    
186
  CASE DEFAULT
187
     Call Die("CheckPeriodicBound"," Iper unknown")
188
  END SELECT
189

    
190

    
191
  if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ==========="
192
END SUBROUTINE CheckPeriodicBound
193

    
194

    
195
 
196

    
197

    
198