← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Util/STM/ol-stm/Src/extrapolate.f

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/STM/ol-stm/Src/extrapolate.f
1
1
! ---
2
 
! Copyright (C) 1996-2016       The SIESTA group
 
2
! Copyright (C) 1996-2016       The SIESTA group
3
3
!  This file is distributed under the terms of the
4
4
!  GNU General Public License: see COPYING in the top directory
5
5
!  or http://www.gnu.org/copyleft/gpl.txt .
46
46
C----------------------------------------------------------
47
47
C INTERNAL VARIABLES
48
48
 
49
 
      INTEGER
50
 
     .  NX, NY, NZ, NX1, NY1, IA1, IA2
 
49
      INTEGER ::  NX, NY, NZ, NX1, NY1, IA1, IA2
51
50
 
52
 
      REAL(DP)
53
 
     .  DECAY, GX, GY, Z, VOLCEL, EAU, V0AU
 
51
      REAL(DP) :: DECAY, GX, GY, Z, VOLCEL, EAU, V0AU
54
52
 
55
53
      REAL(DP), SAVE ::
56
54
     .  STEPZ, A1(3), A2(3), A3(3), VOL, VEC1(3), VEC2(3)
57
55
 
58
 
      COMPLEX(DP)
59
 
     .  EXPSI(0:NPX-1,0:NPY-1)
60
 
 
61
 
      LOGICAL, SAVE :: FIRST
62
 
 
63
 
      EXTERNAL 
64
 
     .  VOLCEL
65
 
 
66
 
!      EXTERNAL
67
 
!     .  dfftw_plan_dft_2d, dfftw_execute, dfftw_destroy_plan
68
 
 
69
 
      DATA FIRST /.TRUE./
70
 
 
 
56
      COMPLEX(DP) :: EXPSI(0:NPX-1,0:NPY-1)
 
57
 
 
58
      LOGICAL, SAVE :: FIRST = .true.
 
59
 
 
60
      EXTERNAL  VOLCEL
71
61
 
72
62
C CHANGE UNITS OF ENERGY TO A.U.
73
63
 
75
65
      V0AU = V0 / 27.21162908D0
76
66
 
77
67
C CHECK THAT STATE IS BOUND (E BELOW VACUUM LEVEL)
78
 
      IF (EAU .LT. V0AU) THEN
 
68
      IF (EAU .GE. V0AU) THEN
 
69
         write(6,*) 'wf NOT considered:'
 
70
         write(6,*) ' ENERGY EIGENVALUE ABOVE VACUUM LEVEL'
 
71
         RETURN
 
72
      endif
79
73
 
80
74
C INITIALIZE VARIABLES ............
81
75
 
82
 
      IF (FIRST) THEN
83
 
 
84
76
! the only problem is the decay
85
77
! which is calculated now by an ABS
86
78
! so in principle this can be commented out
87
 
! N.L. Fevrier 2005
 
79
!     N.L. Fevrier 2005
 
80
 
 
81
!! .... but what is the physical meaning of a decay *towards* the surface?? (AG)
 
82
         
88
83
!       IF (ZMIN .LT. ZREF) THEN
89
84
!         WRITE(6,*) 'error: the reference plane can not be above'
90
85
!         WRITE(6,*) '       the simulation area'
92
87
!         STOP
93
88
!       ENDIF
94
89
 
95
 
C
 
90
      IF (FIRST) THEN
 
91
 
96
92
        IF (NPZ .NE. 1) THEN
97
93
          STEPZ = (ZMAX - ZMIN) / (NPZ - 1)
98
94
        ELSE 
100
96
        ENDIF
101
97
C
102
98
        VOL = VOLCEL(UCELL)
103
 
C
104
 
!!      This was wrong !
105
 
!!        A1=(/UCELL(1,1), UCELL(1,2), UCELL(1,3)/)
106
 
!!        A2=(/UCELL(2,1), UCELL(2,2), UCELL(2,3)/)
107
 
!!        A3=(/UCELL(3,1), UCELL(3,2), UCELL(3,3)/)
108
99
 
109
100
        A1 = UCELL(:,1)
110
101
        A2 = UCELL(:,2)
111
102
        A3 = UCELL(:,3)
112
 
 
113
 
C
114
103
        CALL RECIPROCAL(A1,A2,A3,VOL,VEC1,VEC2)
115
 
C
 
104
 
116
105
        FIRST = .FALSE.
117
 
C
 
106
 
118
107
      ENDIF
119
108
 
120
 
C .....
121
 
 
122
 
C      REWIND(7)
123
 
C      DO NX=0,NPX-1
124
 
C        DO NY=0,NPY-1
125
 
C          WRITE(7,'(4f20.10)') UCELL(1,1)*NX/NPX, 
126
 
C     .             UCELL(2,2)*NY/NPY,
127
 
C     .             DREAL(CW(NX,NY)*DCONJG(CW(NX,NY)))
128
 
C        ENDDO
129
 
C        WRITE(7,*)
130
 
C      ENDDO
131
 
 
132
 
 
133
109
C DO DIRECT FOURIER TRANSFORM TO GET SPATIAL FREQUENCIES OF WF AT 
134
110
C REFERENCE PLANE ........
135
111
 
139
115
      call fftw_execute_dft (plan,cw,cw)
140
116
      call fftw_destroy_plan(plan)
141
117
 
142
 
C .....
143
 
 
144
118
C LOOP OVER SIMULATION HEIGHTS ........
145
119
      plan = fftw_plan_dft_2d (NPY,NPX,EXPSI,EXPSI,FFTW_BACKWARD, 
146
120
     .                        FFTW_ESTIMATE)
 
121
 
147
122
      DO NZ = 1, NPZ
148
123
        Z = ZMIN + (NZ-1)*STEPZ
149
124
 
 
125
        ! NOTE: No "inward propagation...", despite what Nicolas says above
 
126
        ! Points with Z < ZREF are computed from un-propagated wfs
 
127
        IF (Z < ZREF) CYCLE  
 
128
        
150
129
C LOOP OVER POINTS IN XY PLANE ...
151
130
        DO NY = 1, NPY
152
131
          NY1 = NY-1
153
 
C INDEX CHOSEN TO CENTER G's
 
132
C     INDEX CHOSEN TO CENTER G's
 
133
          ! Ambiguous arithmetic here and below (precedence?)
154
134
          IA2 = MOD (NY1,NPY/2)-NY1/(NPY/2)*(NPY/2)
155
135
          DO NX = 1, NPX
156
136
            NX1 = NX-1
160
140
            GY = IA1*VEC1(2) + IA2*VEC2(2)
161
141
 
162
142
C CALCULATE EXPONENT OF DECAY
163
 
           DECAY = DSQRT( 2.0D0*(V0AU-EAU) + 
 
143
           DECAY = SQRT( 2.0D0*(V0AU-EAU) + 
164
144
     .                    (GX+K(1))**2 + (GY+K(2))**2)
165
145
 
166
146
C EXTEND WAVE FUNCTION TO CURRENT HEIGHT
167
 
            EXPSI(NX1,NY1) = CW(NX1,NY1) * DEXP(-DECAY*ABS(Z-ZREF)) / 
 
147
            EXPSI(NX1,NY1) = CW(NX1,NY1) * EXP(-DECAY*ABS(Z-ZREF)) / 
168
148
     .                       (NPX*NPY)
169
149
          ENDDO
170
 
        ENDDO
171
 
C ... END LOOP XY
 
150
        ENDDO        !  LOOP XY
172
151
 
173
152
C DO BACK FOURIER TRANSFORM TO GET REAL SPACE WF AT 
174
153
C REFERENCE PLANE .....
175
154
 
176
155
        call fftw_execute_dft (plan,expsi,expsi)
177
156
 
178
 
C .....
179
 
 
180
157
        CWE(:,:,NZ-1) = EXPSI(:,:)
181
158
 
182
 
      ENDDO
183
 
C ........ END LOOP Z
 
159
      ENDDO     !  LOOP Z
 
160
      
184
161
      call fftw_destroy_plan(plan)
185
162
 
186
 
 
187
 
C      DO NX=0,NPX-1
188
 
C        DO NY=0,NPY-1
189
 
C          WRITE(7,'(3f20.10)') UCELL(1,1)*NX/NPX, 
190
 
C     .             UCELL(2,2)*NY/NPY,
191
 
C     .             DREAL(CW(NX,NY)*DCONJG(CW(NX,NY)))
192
 
C        ENDDO
193
 
C        WRITE(7,*)
194
 
C      ENDDO
195
 
C
196
 
C      REWIND(8)
197
 
C
198
 
C      DO NX=0,NPX-1
199
 
C        DO NY=0,NPY-1
200
 
C          WRITE(8,'(3f20.10)') UCELL(1,1)*NX/NPX, 
201
 
C     .             UCELL(2,2)*NY/NPY,
202
 
C     .             DREAL(CWE(NX,NY,NPZ-1)*DCONJG(CWE(NX,NY,NPZ-1)))
203
 
C        ENDDO
204
 
C        WRITE(8,*)
205
 
C      ENDDO
206
 
 
207
 
 
208
 
 
209
 
        else
210
 
       write(6,*) 'wf NOT considered:'
211
 
       write(6,*) ' ENERGY EIGENVALUE ABOVE VACUUM LEVEL'
212
 
 
213
 
        endif
214
 
 
215
 
        CONTAINS
 
163
      RETURN
 
164
      END
 
165
 
216
166
 
217
167
        SUBROUTINE  reciprocal(a1,a2,a3,vol,vec1,vec2)
218
 
 
219
 
        real(dp)    vol
220
 
        real(dp)    pv(3),a1(3),a2(3),a3(3)
221
 
        real(dp)    vec1(3),vec2(3), PI_D
 
168
        implicit none
 
169
        real*8    vol
 
170
        real*8    pv(3),a1(3),a2(3),a3(3)
 
171
        real*8    vec1(3),vec2(3), PI_D
222
172
 
223
173
        parameter (PI_D = 3.141592653589793238462643383279502884197d0)
224
174
 
230
180
        END SUBROUTINE reciprocal
231
181
 
232
182
        SUBROUTINE pro(v1,v2,pv)
233
 
 
234
 
        real(dp), intent(in)  ::   v1(3),v2(3)
235
 
        real(dp), intent(out) ::   pv(3)
 
183
        implicit none
 
184
        real*8   v1(3),v2(3)
 
185
        real*8   pv(3)
236
186
 
237
187
        pv(1)=v1(2)*v2(3)-v1(3)*v2(2)
238
188
        pv(2)=v1(3)*v2(1)-v1(1)*v2(3)
240
190
 
241
191
        END SUBROUTINE pro
242
192
 
243
 
      END SUBROUTINE EXTRAPOLATE
244
 
 
245
193
 
246
194
 
247
195