1 diff -Nur vasp.4.6.orig/dyna.F vasp.4.6-gamma-fullct/dyna.F
2 --- vasp.4.6.orig/dyna.F 2009-09-16 10:45:06.050116866 +0200
3 +++ vasp.4.6-gamma-fullct/dyna.F 2009-11-05 15:09:26.435422082 +0100
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(NI)/=0._q) THEN
30 + VEC(1)=COS(T_INFO%SD_ROT_Z(NI))*TMP(1)- &
31 + & SIN(T_INFO%SD_ROT_Z(NI))*TMP(2)
32 + VEC(2)=SIN(T_INFO%SD_ROT_Z(NI))*TMP(1)+ &
33 + & COS(T_INFO%SD_ROT_Z(NI))*TMP(2)
35 + IF (T_INFO%SD_ROT_X(NI)/=0._q) THEN
39 + VEC(2)=COS(T_INFO%SD_ROT_X(NI))*TMP(2)- &
40 + & SIN(T_INFO%SD_ROT_X(NI))*TMP(3)
41 + VEC(3)=SIN(T_INFO%SD_ROT_X(NI))*TMP(2)+ &
42 + & COS(T_INFO%SD_ROT_X(NI))*TMP(3)
45 + IF (T_INFO%SD_ROT_X(NI)/=0._q) THEN
49 + VEC(2)=COS(-T_INFO%SD_ROT_X(NI))*TMP(2)- &
50 + & SIN(-T_INFO%SD_ROT_X(NI))*TMP(3)
51 + VEC(3)=SIN(-T_INFO%SD_ROT_X(NI))*TMP(2)+ &
52 + & COS(-T_INFO%SD_ROT_X(NI))*TMP(3)
54 + IF (T_INFO%SD_ROT_Z(NI)/=0._q) THEN
58 + VEC(1)=COS(-T_INFO%SD_ROT_Z(NI))*TMP(1)- &
59 + & SIN(-T_INFO%SD_ROT_Z(NI))*TMP(2)
60 + VEC(2)=SIN(-T_INFO%SD_ROT_Z(NI))*TMP(1)+ &
61 + & COS(-T_INFO%SD_ROT_Z(NI))*TMP(2)
71 +!*********************************************************************
73 ! subroutine to set selected force components to zero
75 !*********************************************************************
76 @@ -1184,19 +1249,44 @@
77 IF (T_INFO%LSDYN) THEN
78 !-----Selective dynamics here ... :
80 +!-----transformed selective dynamics
81 + IF(T_INFO%SD_ROT) THEN
85 + CALL SD_ROT(T_INFO,VTMP,'F',NI)
90 !-----reset the velocities of the selected coordinates ...
92 IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
94 +!-----transform back to direct coordinates
95 + IF(T_INFO%SD_ROT) THEN
99 + CALL SD_ROT(T_INFO,VTMP,'B',NI)
104 !-----and reset the selected force coordinates (warning: forces are
105 ! given in cartesian, selection is made in direct coordinates!):
109 CALL KARDIR(1,VTMP,LATT_CUR%B)
110 +!-----transformed selective dynamics
111 + IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'F',NI)
112 +!-----the actual reset
114 IF (.NOT.T_INFO%LSFOR(M,NI)) VTMP(M)=0._q
116 +!-----transform back to direct coordinates
117 + IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'B',NI)
118 CALL DIRKAR(1,VTMP,LATT_CUR%A)
121 @@ -1222,10 +1312,30 @@
122 IF (T_INFO%LSDYN) THEN
123 !-----Selective dynamics here ... :
125 +!-----transformed selective dynamics
126 + IF(T_INFO%SD_ROT) THEN
130 + CALL SD_ROT(T_INFO,VTMP,'F',NI)
135 !-----reset the velocities of the selected coordinates ...
137 IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
139 +!-----transform back to direct coordinates
140 + IF(T_INFO%SD_ROT) THEN
144 + CALL SD_ROT(T_INFO,VTMP,'B',NI)
152 diff -Nur vasp.4.6.orig/main.F vasp.4.6-gamma-fullct/main.F
153 --- vasp.4.6.orig/main.F 2009-09-16 10:45:05.747597375 +0200
154 +++ vasp.4.6-gamma-fullct/main.F 2009-11-05 15:09:26.449951618 +0100
155 @@ -1929,7 +1929,8 @@
157 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
158 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
159 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
160 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
161 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
163 ! write AE charge density
164 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
165 @@ -1951,7 +1952,8 @@
167 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
168 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
169 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
170 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
171 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
173 ! write AE core charge density
174 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
175 @@ -3478,7 +3480,7 @@
178 CALL OUTPOS(13,.TRUE.,T_INFO%SZNAM2,LATT_CUR%SCALE,DYN%AC,T_INFO%NTYP,T_INFO%NITYP,T_INFO%LSDYN, &
179 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR )
180 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
181 CALL OUTPOS_TRAIL(13,IO%LOPEN, LATT_CUR, T_INFO, DYN)
184 @@ -3489,7 +3491,7 @@
187 CALL OUTPOS(70,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
188 - & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR)
189 + & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
193 @@ -3828,7 +3830,7 @@
196 CALL OUTPOS(18,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
197 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
198 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
200 ! if you uncomment the following lines the pseudo core charge density
201 ! is added to the pseudo charge density
202 @@ -3855,7 +3857,7 @@
203 IF (IO%LOPEN) OPEN(IO%IUVTOT,FILE='LOCPOT',STATUS='UNKNOWN')
205 CALL OUTPOS(IO%IUVTOT,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
206 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
207 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
209 ! comment out the following line to add exchange correlation
211 @@ -3936,7 +3938,8 @@
213 CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
214 & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
215 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
216 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
217 + & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
219 ! write AE charge density
220 CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
221 @@ -3959,7 +3962,7 @@
223 OPEN(UNIT=53,FILE='ELFCAR',STATUS='UNKNOWN')
224 CALL OUTPOS(53,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
225 - & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
226 + & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
230 diff -Nur vasp.4.6.orig/pardens.F vasp.4.6-gamma-fullct/pardens.F
231 --- vasp.4.6.orig/pardens.F 2009-09-16 10:45:06.294577287 +0200
232 +++ vasp.4.6-gamma-fullct/pardens.F 2009-11-05 15:09:26.453773187 +0100
234 ! I do not call <<write_pard>> here, because I handle the total charge a little bit
235 ! differnt, and the output to <<CHG>> is also differnt
236 do_io CALL OUTPOS(iuchgcar,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP, &
237 - T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
238 + T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
240 CALL OUTCHG(GRIDC,iuchgcar,.TRUE.,CHTOT)
243 do_io CLOSE(iuchgcar)
245 do_io CALL OUTPOS(iuchg,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP, &
246 - T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
247 + T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
249 CALL OUTCHG(GRIDC,iuchg,.FALSE.,CHTOT)
251 @@ -1154,7 +1154,7 @@
253 ! and use the usual output routines
254 CALL OUTPOS(iunit,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP, &
255 - T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
256 + T_INFO%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
259 CALL OUTCHG(GRIDC,iunit,.FALSE.,CHTOT)
260 diff -Nur vasp.4.6.orig/poscar.F vasp.4.6-gamma-fullct/poscar.F
261 --- vasp.4.6.orig/poscar.F 2009-09-16 10:45:06.074097686 +0200
262 +++ vasp.4.6-gamma-fullct/poscar.F 2009-11-05 15:19:56.505775347 +0100
264 INTEGER I,NT,NI,NSCALE
265 REAL(q) SCALEX,SCALEY,SCALEZ
267 + REAL(q),PARAMETER :: PI=3.1415926536_q
269 OPEN(UNIT=15,FILE=DIR_APP(1:DIR_LEN)//'POSCAR',STATUS='OLD')
275 + T_INFO%SD_ROT=.FALSE.
276 + IF ((CSEL=='T').OR.(CSEL=='t').OR. &
277 + & (CSEL=='R').OR.(CSEL=='r')) THEN
278 +!-----transformed/rotated selective dynamics switched on
279 + T_INFO%SD_ROT=.TRUE.
280 +!-----print some info and convert to rad
281 + WRITE(IU6,*) ' rotated selective dynamics switched on'
282 +!-----switch on selective dynamics for further processing
286 T_INFO%LSDYN=((CSEL=='S').OR.(CSEL=='s'))
287 IF (T_INFO%LSDYN) READ(15,'(A1)') CSEL
288 IF (CSEL=='K'.OR.CSEL=='k'.OR. &
290 & DYN%D2C(3,NIOND), &
291 & DYN%VEL(3,NIOND),DYN%D2(3,NIOND),DYN%D3(3,NIOND))
293 + ALLOCATE(T_INFO%SD_ROT_Z(NIOND))
294 + ALLOCATE(T_INFO%SD_ROT_X(NIOND))
296 + T_INFO%SD_ROT_Z=0._q
297 + T_INFO%SD_ROT_X=0._q
299 ! alias T_INFO%POSION
300 T_INFO%POSION => DYN%POSION
302 @@ -298,12 +316,21 @@
306 + IF (T_INFO%SD_ROT) THEN
307 + READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI), &
308 + & T_INFO%LSFOR(1,NI),T_INFO%LSFOR(2,NI),T_INFO%LSFOR(3,NI), &
309 + & T_INFO%SD_ROT_Z(NI), T_INFO%SD_ROT_X(NI)
310 +!-----convert to rad
311 + T_INFO%SD_ROT_Z(NI)=PI*T_INFO%SD_ROT_Z(NI)/180.0
312 + T_INFO%SD_ROT_X(NI)=PI*T_INFO%SD_ROT_X(NI)/180.0
314 IF (T_INFO%LSDYN) THEN
315 READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI), &
316 & T_INFO%LSFOR(1,NI),T_INFO%LSFOR(2,NI),T_INFO%LSFOR(3,NI)
318 READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI)
324 @@ -538,17 +565,19 @@
325 !***********************************************************************
327 SUBROUTINE OUTPOS(IU,LLONG,SZNAM,SCALE,A,NTYP,NITYP,LSDYN, &
328 - & NIONS,POSION,LSFOR )
329 + & NIONS,POSION,LSFOR,SD_ROT,SD_ROT_Z,SD_ROT_X)
331 IMPLICIT REAL(q) (A-H,O-Z)
333 - LOGICAL LLONG,LSDYN,LSFOR
334 + LOGICAL LLONG,LSDYN,LSFOR,SD_ROT
335 + REAL(q), POINTER :: SD_ROT_Z(:), SD_ROT_X(:)
337 DIMENSION POSION(3,NIONS)
338 DIMENSION LSFOR (3,NIONS)
339 DIMENSION NITYP(NTYP)
342 + REAL(q),PARAMETER :: PI=3.1415926536_q
345 WRITE(IU,'(A40)') SZNAM
346 @@ -562,14 +591,26 @@
347 WRITE(IU,FORM) (A(1,I)/SCALE,A(2,I)/SCALE,A(3,I)/SCALE,I=1,3)
349 WRITE(IU,'(20I4)') (NITYP(NT),NT=1,NTYP)
350 - IF (LSDYN) WRITE(13,'(A18)') 'Selective dynamics'
353 + WRITE(13,'(A30)') 'Transformed selective dynamics'
355 + WRITE(13,'(A18)') 'Selective dynamics'
358 WRITE(IU,'(A6)')'Direct'
364 + FORM='(3F20.16,3L4,2F20.12)'
369 + FORM='(3F10.6,3L2,2F10.2)'
380 + & (POSION(1,NI),POSION(2,NI),POSION(3,NI), &
381 + & LSFOR(1,NI),LSFOR(2,NI),LSFOR(3,NI), &
382 + & SD_ROT_Z(NI)*180.0/PI, SD_ROT_X(NI)*180.0/PI,NI=1,NIONS)
386 & (POSION(1,NI),POSION(2,NI),POSION(3,NI), &
389 & (POSION(1,NI),POSION(2,NI),POSION(3,NI),NI=1,NIONS)
392 IF (.NOT.LLONG) WRITE(IU,*)
395 diff -Nur vasp.4.6.orig/poscar.inc vasp.4.6-gamma-fullct/poscar.inc
396 --- vasp.4.6.orig/poscar.inc 2009-09-16 10:45:06.791087379 +0200
397 +++ vasp.4.6-gamma-fullct/poscar.inc 2009-11-03 11:58:39.388412198 +0100
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), POINTER :: SD_ROT_Z(:) ! constraints rotation (z-axis, 1st rot)
404 + REAL(q), POINTER :: 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