1 diff -Nur vasp.4.6/dyna.F vasp.4.6-ct/dyna.F
2 --- vasp.4.6/dyna.F 2009-07-22 15:12:20.032980043 +0200
3 +++ vasp.4.6-ct/dyna.F 2009-09-18 08:32:08.666016928 +0200
6 !*********************************************************************
8 +! subroutine to transform vectors into the rotated constraints system
10 +!*********************************************************************
12 + SUBROUTINE SD_ROT(T_INFO,VEC,DIR,NI)
16 + TYPE (type_info) T_INFO
17 + REAL(q) VEC(3),TMP(3)
22 + IF((.NOT. T_INFO%LSFOR(1,NI)).OR. &
23 + & (.NOT. T_INFO%LSFOR(2,NI)).OR. &
24 + & (.NOT. T_INFO%LSFOR(3,NI))) THEN
26 + IF (T_INFO%SD_ROT_Z/=0._q) THEN
30 + VEC(1)=COS(T_INFO%SD_ROT_Z)*TMP(1)- &
31 + & SIN(T_INFO%SD_ROT_Z)*TMP(2)
32 + VEC(2)=SIN(T_INFO%SD_ROT_Z)*TMP(1)+ &
33 + & COS(T_INFO%SD_ROT_Z)*TMP(2)
34 + WRITE(*,*) 'SD rotate DEBUG (Z^-1): ',NI
35 + WRITE(*,*) ' F: ', TMP(1),' ',TMP(2),' ',TMP(3)
36 + WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
38 + IF (T_INFO%SD_ROT_X/=0._q) THEN
42 + VEC(2)=COS(T_INFO%SD_ROT_X)*TMP(2)- &
43 + & SIN(T_INFO%SD_ROT_X)*TMP(3)
44 + VEC(3)=SIN(T_INFO%SD_ROT_X)*TMP(2)+ &
45 + & COS(T_INFO%SD_ROT_X)*TMP(3)
46 + WRITE(*,*) 'SD rotate DEBUG (X^-1): ',NI
47 + WRITE(*,*) ' F: ', TMP(1),' ',TMP(2),' ',TMP(3)
48 + WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
51 + IF (T_INFO%SD_ROT_X/=0._q) THEN
55 + VEC(2)=COS(-T_INFO%SD_ROT_X)*TMP(2)- &
56 + & SIN(-T_INFO%SD_ROT_X)*TMP(3)
57 + VEC(3)=SIN(-T_INFO%SD_ROT_X)*TMP(2)+ &
58 + & COS(-T_INFO%SD_ROT_X)*TMP(3)
59 + WRITE(*,*) 'SD rotate DEBUG (X): ',NI
60 + WRITE(*,*) ' B: ', TMP(1),' ',TMP(2),' ',TMP(3)
61 + WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
63 + IF (T_INFO%SD_ROT_Z/=0._q) THEN
67 + VEC(1)=COS(-T_INFO%SD_ROT_Z)*TMP(1)- &
68 + & SIN(-T_INFO%SD_ROT_Z)*TMP(2)
69 + VEC(2)=SIN(-T_INFO%SD_ROT_Z)*TMP(1)+ &
70 + & COS(-T_INFO%SD_ROT_Z)*TMP(2)
71 + WRITE(*,*) 'SD rotate DEBUG (Z): ',NI
72 + WRITE(*,*) ' B: ', TMP(1),' ',TMP(2),' ',TMP(3)
73 + WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
83 +!*********************************************************************
85 ! subroutine to set selected force components to zero
87 !*********************************************************************
88 @@ -1184,19 +1261,44 @@
89 IF (T_INFO%LSDYN) THEN
90 !-----Selective dynamics here ... :
92 +!-----transformed selective dynamics
93 + IF(T_INFO%SD_ROT) THEN
97 + CALL SD_ROT(T_INFO,VTMP,'F',NI)
102 !-----reset the velocities of the selected coordinates ...
104 IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
106 +!-----transform back to direct coordinates
107 + IF(T_INFO%SD_ROT) THEN
111 + CALL SD_ROT(T_INFO,VTMP,'B',NI)
116 !-----and reset the selected force coordinates (warning: forces are
117 ! given in cartesian, selection is made in direct coordinates!):
121 CALL KARDIR(1,VTMP,LATT_CUR%B)
122 +!-----transformed selective dynamics
123 + IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'F',NI)
124 +!-----the actual reset
126 IF (.NOT.T_INFO%LSFOR(M,NI)) VTMP(M)=0._q
128 +!-----transform back to direct coordinates
129 + IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'B',NI)
130 CALL DIRKAR(1,VTMP,LATT_CUR%A)
133 @@ -1222,10 +1324,30 @@
134 IF (T_INFO%LSDYN) THEN
135 !-----Selective dynamics here ... :
137 +!-----transformed selective dynamics
138 + IF(T_INFO%SD_ROT) THEN
142 + CALL SD_ROT(T_INFO,VTMP,'F',NI)
147 !-----reset the velocities of the selected coordinates ...
149 IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
151 +!-----transform back to direct coordinates
152 + IF(T_INFO%SD_ROT) THEN
156 + CALL SD_ROT(T_INFO,VTMP,'B',NI)
164 diff -Nur vasp.4.6/main.F vasp.4.6-ct/main.F
165 --- vasp.4.6/main.F 2009-07-22 15:12:20.044978904 +0200
166 +++ vasp.4.6-ct/main.F 2009-09-18 08:55:20.595465366 +0200
167 @@ -1929,7 +1929,8 @@
169 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
170 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
171 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
172 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
173 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
175 ! write AE charge density
176 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
177 @@ -1951,7 +1952,8 @@
179 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
180 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
181 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
182 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
183 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
185 ! write AE core charge density
186 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
187 @@ -3478,7 +3480,7 @@
190 CALL OUTPOS(13,.TRUE.,T_INFO%SZNAM2,LATT_CUR%SCALE,DYN%AC,T_INFO%NTYP,T_INFO%NITYP,T_INFO%LSDYN, &
191 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR )
192 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
193 CALL OUTPOS_TRAIL(13,IO%LOPEN, LATT_CUR, T_INFO, DYN)
196 @@ -3489,7 +3491,7 @@
199 CALL OUTPOS(70,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
200 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR)
201 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
205 @@ -3828,7 +3830,7 @@
208 CALL OUTPOS(18,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
209 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
210 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
212 ! if you uncomment the following lines the pseudo core charge density
213 ! is added to the pseudo charge density
214 @@ -3855,7 +3857,7 @@
215 IF (IO%LOPEN) OPEN(IO%IUVTOT,FILE='LOCPOT',STATUS='UNKNOWN')
217 CALL OUTPOS(IO%IUVTOT,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
218 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
219 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
221 ! comment out the following line to add exchange correlation
223 @@ -3936,7 +3938,8 @@
225 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
226 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
227 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
228 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
229 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
231 ! write AE charge density
232 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
233 @@ -3959,7 +3962,7 @@
235 OPEN(UNIT=53,FILE='ELFCAR',STATUS='UNKNOWN')
236 CALL OUTPOS(53,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
237 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
238 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
242 diff -Nur vasp.4.6/main.f90 vasp.4.6-ct/main.f90
243 --- vasp.4.6/main.f90 2009-07-22 15:14:38.603841627 +0200
244 +++ vasp.4.6-ct/main.f90 2009-09-18 08:55:22.638156662 +0200
245 @@ -1830,7 +1830,8 @@
247 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
248 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
249 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
250 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
251 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
253 ! write AE charge density
254 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
255 @@ -1852,7 +1853,8 @@
257 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
258 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
259 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
260 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
261 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
263 ! write AE core charge density
264 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
265 @@ -3369,7 +3371,7 @@
268 CALL OUTPOS(13,.TRUE.,T_INFO%SZNAM2,LATT_CUR%SCALE,DYN%AC,T_INFO%NTYP,T_INFO%NITYP,T_INFO%LSDYN, &
269 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR )
270 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
271 CALL OUTPOS_TRAIL(13,IO%LOPEN, LATT_CUR, T_INFO, DYN)
274 @@ -3380,7 +3382,7 @@
277 CALL OUTPOS(70,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
278 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR)
279 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
283 @@ -3712,7 +3714,7 @@
286 CALL OUTPOS(18,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
287 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
288 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
290 ! if you uncomment the following lines the pseudo core charge density
291 ! is added to the pseudo charge density
292 @@ -3739,7 +3741,7 @@
293 IF (IO%LOPEN) OPEN(IO%IUVTOT,FILE='LOCPOT',STATUS='UNKNOWN')
295 CALL OUTPOS(IO%IUVTOT,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
296 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
297 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
299 ! comment out the following line to add exchange correlation
301 @@ -3820,7 +3822,8 @@
303 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
304 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
305 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
306 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
307 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
309 ! write AE charge density
310 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
311 @@ -3843,7 +3846,7 @@
313 OPEN(UNIT=53,FILE='ELFCAR',STATUS='UNKNOWN')
314 CALL OUTPOS(53,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
315 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
316 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
320 diff -Nur vasp.4.6/poscar.F vasp.4.6-ct/poscar.F
321 --- vasp.4.6/poscar.F 2009-07-22 15:12:20.052978145 +0200
322 +++ vasp.4.6-ct/poscar.F 2009-09-18 09:17:37.370283715 +0200
324 INTEGER I,NT,NI,NSCALE
325 REAL(q) SCALEX,SCALEY,SCALEZ
327 + REAL(q),PARAMETER :: PI=3.1415926536_q
329 OPEN(UNIT=15,FILE=DIR_APP(1:DIR_LEN)//'POSCAR',STATUS='OLD')
335 + T_INFO%SD_ROT=.FALSE.
336 + IF ((CSEL=='T').OR.(CSEL=='t').OR. &
337 + & (CSEL=='R').OR.(CSEL=='r')) THEN
338 +!-----transformed/rotated selective dynamics switched on
339 +!-----read in rotation angles
340 + T_INFO%SD_ROT=.TRUE.
341 + T_INFO%SD_ROT_Z=0._q
342 + T_INFO%SD_ROT_X=0._q
343 + READ(15,*,ERR=400,END=400) T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X
344 +!-----print some info and convert to rad
345 + WRITE(IU6,*) ' rotated selective dynamics:'
346 + WRITE(IU6,*) ' Z: ', T_INFO%SD_ROT_Z
347 + WRITE(IU6,*) ' rX: ', T_INFO%SD_ROT_X
348 + T_INFO%SD_ROT_Z=PI*T_INFO%SD_ROT_Z/180.0
349 + T_INFO%SD_ROT_X=PI*T_INFO%SD_ROT_X/180.0
350 +!-----switch on selective dynamics for further processing
354 T_INFO%LSDYN=((CSEL=='S').OR.(CSEL=='s'))
355 IF (T_INFO%LSDYN) READ(15,'(A1)') CSEL
356 IF (CSEL=='K'.OR.CSEL=='k'.OR. &
357 @@ -538,17 +558,19 @@
358 !***********************************************************************
360 SUBROUTINE OUTPOS(IU,LLONG,SZNAM,SCALE,A,NTYP,NITYP,LSDYN, &
361 - & NIONS,POSION,LSFOR )
362 + & NIONS,POSION,LSFOR,SD_ROT,SD_ROT_Z,SD_ROT_X)
364 IMPLICIT REAL(q) (A-H,O-Z)
366 - LOGICAL LLONG,LSDYN,LSFOR
367 + LOGICAL LLONG,LSDYN,LSFOR,SD_ROT
368 + REAL(q) SD_ROT_Z,SD_ROT_X
370 DIMENSION POSION(3,NIONS)
371 DIMENSION LSFOR (3,NIONS)
372 DIMENSION NITYP(NTYP)
375 + REAL(q),PARAMETER :: PI=3.1415926536_q
378 WRITE(IU,'(A40)') SZNAM
380 WRITE(IU,FORM) (A(1,I)/SCALE,A(2,I)/SCALE,A(3,I)/SCALE,I=1,3)
382 WRITE(IU,'(20I4)') (NITYP(NT),NT=1,NTYP)
383 - IF (LSDYN) WRITE(13,'(A18)') 'Selective dynamics'
386 + WRITE(13,'(A30)') 'Transformed selective dynamics'
387 + WRITE(13,'(2F)') SD_ROT_Z*180.0/PI,SD_ROT_X*180.0/PI
389 + WRITE(13,'(A18)') 'Selective dynamics'
392 WRITE(IU,'(A6)')'Direct'
395 diff -Nur vasp.4.6/poscar.inc vasp.4.6-ct/poscar.inc
396 --- vasp.4.6/poscar.inc 2009-07-22 15:05:24.076290013 +0200
397 +++ vasp.4.6-ct/poscar.inc 2009-09-18 08:01:29.234012116 +0200
399 INTEGER NIONS ! actual number of ions
400 INTEGER NIONP ! actual number of ions inc. empty spheres
401 LOGICAL LSDYN ! selective dynamics (yes/ no)
402 + LOGICAL SD_ROT ! rotated selective dynamics (yes/ no)
403 + REAL(q) SD_ROT_Z ! constraints rotation (z-axis, 1st rot)
404 + REAL(q) SD_ROT_X ! constraints rotation (x-axis, 2nd rot)
405 LOGICAL LDIRCO ! positions in direct/recproc. lattice
406 REAL(q), POINTER :: POSION(:,:) ! positions usually same as DYN%POSION
407 LOGICAL,POINTER :: LSFOR(:,:) ! selective dynamics