--- /dev/null
+diff -Nur vasp.4.6/dyna.F vasp.4.6-ct/dyna.F
+--- vasp.4.6/dyna.F 2009-07-22 15:12:20.032980043 +0200
++++ vasp.4.6-ct/dyna.F 2009-09-18 08:32:08.666016928 +0200
+@@ -1164,6 +1164,83 @@
+
+ !*********************************************************************
+ !
++! subroutine to transform vectors into the rotated constraints system
++!
++!*********************************************************************
++
++ SUBROUTINE SD_ROT(T_INFO,VEC,DIR,NI)
++ USE poscar
++ IMPLICIT NONE
++
++ TYPE (type_info) T_INFO
++ REAL(q) VEC(3),TMP(3)
++ CHARACTER*1 DIR
++
++ INTEGER M,NI
++
++ IF((.NOT. T_INFO%LSFOR(1,NI)).OR. &
++ & (.NOT. T_INFO%LSFOR(2,NI)).OR. &
++ & (.NOT. T_INFO%LSFOR(3,NI))) THEN
++ IF (DIR=='F') THEN
++ IF (T_INFO%SD_ROT_Z/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(1)=COS(T_INFO%SD_ROT_Z)*TMP(1)- &
++ & SIN(T_INFO%SD_ROT_Z)*TMP(2)
++ VEC(2)=SIN(T_INFO%SD_ROT_Z)*TMP(1)+ &
++ & COS(T_INFO%SD_ROT_Z)*TMP(2)
++ WRITE(*,*) 'SD rotate DEBUG (Z^-1): ',NI
++ WRITE(*,*) ' F: ', TMP(1),' ',TMP(2),' ',TMP(3)
++ WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
++ ENDIF
++ IF (T_INFO%SD_ROT_X/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(2)=COS(T_INFO%SD_ROT_X)*TMP(2)- &
++ & SIN(T_INFO%SD_ROT_X)*TMP(3)
++ VEC(3)=SIN(T_INFO%SD_ROT_X)*TMP(2)+ &
++ & COS(T_INFO%SD_ROT_X)*TMP(3)
++ WRITE(*,*) 'SD rotate DEBUG (X^-1): ',NI
++ WRITE(*,*) ' F: ', TMP(1),' ',TMP(2),' ',TMP(3)
++ WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
++ ENDIF
++ ELSE
++ IF (T_INFO%SD_ROT_X/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(2)=COS(-T_INFO%SD_ROT_X)*TMP(2)- &
++ & SIN(-T_INFO%SD_ROT_X)*TMP(3)
++ VEC(3)=SIN(-T_INFO%SD_ROT_X)*TMP(2)+ &
++ & COS(-T_INFO%SD_ROT_X)*TMP(3)
++ WRITE(*,*) 'SD rotate DEBUG (X): ',NI
++ WRITE(*,*) ' B: ', TMP(1),' ',TMP(2),' ',TMP(3)
++ WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
++ ENDIF
++ IF (T_INFO%SD_ROT_Z/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(1)=COS(-T_INFO%SD_ROT_Z)*TMP(1)- &
++ & SIN(-T_INFO%SD_ROT_Z)*TMP(2)
++ VEC(2)=SIN(-T_INFO%SD_ROT_Z)*TMP(1)+ &
++ & COS(-T_INFO%SD_ROT_Z)*TMP(2)
++ WRITE(*,*) 'SD rotate DEBUG (Z): ',NI
++ WRITE(*,*) ' B: ', TMP(1),' ',TMP(2),' ',TMP(3)
++ WRITE(*,*) ' to: ',VEC(1),' ',VEC(2),' ',VEC(3)
++ ENDIF
++ ENDIF
++ ENDIF
++
++ RETURN
++
++ END
++
++
++!*********************************************************************
++!
+ ! subroutine to set selected force components to zero
+ !
+ !*********************************************************************
+@@ -1184,19 +1261,44 @@
+ IF (T_INFO%LSDYN) THEN
+ !-----Selective dynamics here ... :
+ DO NI=1,T_INFO%NIONS
++!-----transformed selective dynamics
++ IF(T_INFO%SD_ROT) THEN
++ DO M=1,3
++ VTMP(M)=VEL(M,NI)
++ ENDDO
++ CALL SD_ROT(T_INFO,VTMP,'F',NI)
++ DO M=1,3
++ VEL(M,NI)=VTMP(M)
++ ENDDO
++ ENDIF
+ !-----reset the velocities of the selected coordinates ...
+ DO M=1,3
+ IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
+ ENDDO
++!-----transform back to direct coordinates
++ IF(T_INFO%SD_ROT) THEN
++ DO M=1,3
++ VTMP(M)=VEL(M,NI)
++ ENDDO
++ CALL SD_ROT(T_INFO,VTMP,'B',NI)
++ DO M=1,3
++ VEL(M,NI)=VTMP(M)
++ ENDDO
++ ENDIF
+ !-----and reset the selected force coordinates (warning: forces are
+ ! given in cartesian, selection is made in direct coordinates!):
+ DO M=1,3
+ VTMP(M)=TIFOR(M,NI)
+ ENDDO
+ CALL KARDIR(1,VTMP,LATT_CUR%B)
++!-----transformed selective dynamics
++ IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'F',NI)
++!-----the actual reset
+ DO M=1,3
+ IF (.NOT.T_INFO%LSFOR(M,NI)) VTMP(M)=0._q
+ ENDDO
++!-----transform back to direct coordinates
++ IF(T_INFO%SD_ROT) CALL SD_ROT(T_INFO,VTMP,'B',NI)
+ CALL DIRKAR(1,VTMP,LATT_CUR%A)
+ DO M=1,3
+ TIFOR(M,NI)=VTMP(M)
+@@ -1222,10 +1324,30 @@
+ IF (T_INFO%LSDYN) THEN
+ !-----Selective dynamics here ... :
+ DO NI=1,T_INFO%NIONS
++!-----transformed selective dynamics
++ IF(T_INFO%SD_ROT) THEN
++ DO M=1,3
++ VTMP(M)=VEL(M,NI)
++ ENDDO
++ CALL SD_ROT(T_INFO,VTMP,'F',NI)
++ DO M=1,3
++ VEL(M,NI)=VTMP(M)
++ ENDDO
++ ENDIF
+ !-----reset the velocities of the selected coordinates ...
+ DO M=1,3
+ IF (.NOT.T_INFO%LSFOR(M,NI)) VEL(M,NI)=0._q
+ ENDDO
++!-----transform back to direct coordinates
++ IF(T_INFO%SD_ROT) THEN
++ DO M=1,3
++ VTMP(M)=VEL(M,NI)
++ ENDDO
++ CALL SD_ROT(T_INFO,VTMP,'B',NI)
++ DO M=1,3
++ VEL(M,NI)=VTMP(M)
++ ENDDO
++ ENDIF
+ ENDDO
+ ENDIF
+ END SUBROUTINE
+diff -Nur vasp.4.6/main.F vasp.4.6-ct/main.F
+--- vasp.4.6/main.F 2009-07-22 15:12:20.044978904 +0200
++++ vasp.4.6-ct/main.F 2009-09-18 08:55:20.595465366 +0200
+@@ -1929,7 +1929,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+ ! write AE charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -1951,7 +1952,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+ ! write AE core charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -3478,7 +3480,7 @@
+ io_begin
+
+ CALL OUTPOS(13,.TRUE.,T_INFO%SZNAM2,LATT_CUR%SCALE,DYN%AC,T_INFO%NTYP,T_INFO%NITYP,T_INFO%LSDYN, &
+- & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR )
++ & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ CALL OUTPOS_TRAIL(13,IO%LOPEN, LATT_CUR, T_INFO, DYN)
+
+ io_end
+@@ -3489,7 +3491,7 @@
+
+ io_begin
+ CALL OUTPOS(70,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+
+ DO ISP=1,WDES%NCDIJ
+@@ -3828,7 +3830,7 @@
+ REWIND 18
+ ! since t
+ CALL OUTPOS(18,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+ ! if you uncomment the following lines the pseudo core charge density
+ ! is added to the pseudo charge density
+@@ -3855,7 +3857,7 @@
+ IF (IO%LOPEN) OPEN(IO%IUVTOT,FILE='LOCPOT',STATUS='UNKNOWN')
+ REWIND IO%IUVTOT
+ CALL OUTPOS(IO%IUVTOT,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+ ! comment out the following line to add exchange correlation
+ INFO%LEXCHG=-1
+@@ -3936,7 +3938,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+ ! write AE charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -3959,7 +3962,7 @@
+ io_begin
+ OPEN(UNIT=53,FILE='ELFCAR',STATUS='UNKNOWN')
+ CALL OUTPOS(53,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+
+ DO ISP=1,WDES%NCDIJ
+diff -Nur vasp.4.6/main.f90 vasp.4.6-ct/main.f90
+--- vasp.4.6/main.f90 2009-07-22 15:14:38.603841627 +0200
++++ vasp.4.6-ct/main.f90 2009-09-18 08:55:22.638156662 +0200
+@@ -1830,7 +1830,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ ! write AE charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -1852,7 +1853,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ ! write AE core charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -3369,7 +3371,7 @@
+
+
+ CALL OUTPOS(13,.TRUE.,T_INFO%SZNAM2,LATT_CUR%SCALE,DYN%AC,T_INFO%NTYP,T_INFO%NITYP,T_INFO%LSDYN, &
+- & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR )
++ & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ CALL OUTPOS_TRAIL(13,IO%LOPEN, LATT_CUR, T_INFO, DYN)
+
+
+@@ -3380,7 +3382,7 @@
+
+
+ CALL OUTPOS(70,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSIOC,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+
+ DO ISP=1,WDES%NCDIJ
+@@ -3712,7 +3714,7 @@
+ REWIND 18
+ ! since t
+ CALL OUTPOS(18,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ ! if you uncomment the following lines the pseudo core charge density
+ ! is added to the pseudo charge density
+@@ -3739,7 +3741,7 @@
+ IF (IO%LOPEN) OPEN(IO%IUVTOT,FILE='LOCPOT',STATUS='UNKNOWN')
+ REWIND IO%IUVTOT
+ CALL OUTPOS(IO%IUVTOT,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ ! comment out the following line to add exchange correlation
+ INFO%LEXCHG=-1
+@@ -3820,7 +3822,8 @@
+ ! write header
+ CALL OUTPOS(99,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A, &
+ & T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR, &
++ & T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ ! write AE charge density
+ CALL OUTCHG(GRIDC,99,.TRUE.,CHTOT)
+@@ -3843,7 +3846,7 @@
+
+ OPEN(UNIT=53,FILE='ELFCAR',STATUS='UNKNOWN')
+ CALL OUTPOS(53,.FALSE.,INFO%SZNAM1,LATT_CUR%SCALE,LATT_CUR%A,T_INFO%NTYP,T_INFO%NITYP,.FALSE., &
+- & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR)
++ & T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+
+ DO ISP=1,WDES%NCDIJ
+diff -Nur vasp.4.6/poscar.F vasp.4.6-ct/poscar.F
+--- vasp.4.6/poscar.F 2009-07-22 15:12:20.052978145 +0200
++++ vasp.4.6-ct/poscar.F 2009-09-18 09:17:37.370283715 +0200
+@@ -181,6 +181,7 @@
+ INTEGER I,NT,NI,NSCALE
+ REAL(q) SCALEX,SCALEY,SCALEZ
+ REAL(q) POTIMR
++ REAL(q),PARAMETER :: PI=3.1415926536_q
+
+ OPEN(UNIT=15,FILE=DIR_APP(1:DIR_LEN)//'POSCAR',STATUS='OLD')
+
+@@ -269,6 +270,25 @@
+ ENDIF
+
+ READ(15,'(A1)') CSEL
++ T_INFO%SD_ROT=.FALSE.
++ IF ((CSEL=='T').OR.(CSEL=='t').OR. &
++ & (CSEL=='R').OR.(CSEL=='r')) THEN
++!-----transformed/rotated selective dynamics switched on
++!-----read in rotation angles
++ T_INFO%SD_ROT=.TRUE.
++ T_INFO%SD_ROT_Z=0._q
++ T_INFO%SD_ROT_X=0._q
++ READ(15,*,ERR=400,END=400) T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X
++!-----print some info and convert to rad
++ WRITE(IU6,*) ' rotated selective dynamics:'
++ WRITE(IU6,*) ' Z: ', T_INFO%SD_ROT_Z
++ WRITE(IU6,*) ' rX: ', T_INFO%SD_ROT_X
++ T_INFO%SD_ROT_Z=PI*T_INFO%SD_ROT_Z/180.0
++ T_INFO%SD_ROT_X=PI*T_INFO%SD_ROT_X/180.0
++!-----switch on selective dynamics for further processing
++ CSEL='s'
++ ENDIF
++
+ T_INFO%LSDYN=((CSEL=='S').OR.(CSEL=='s'))
+ IF (T_INFO%LSDYN) READ(15,'(A1)') CSEL
+ IF (CSEL=='K'.OR.CSEL=='k'.OR. &
+@@ -538,17 +558,19 @@
+ !***********************************************************************
+
+ SUBROUTINE OUTPOS(IU,LLONG,SZNAM,SCALE,A,NTYP,NITYP,LSDYN, &
+- & NIONS,POSION,LSFOR )
++ & NIONS,POSION,LSFOR,SD_ROT,SD_ROT_Z,SD_ROT_X)
+ USE prec
+ IMPLICIT REAL(q) (A-H,O-Z)
+
+- LOGICAL LLONG,LSDYN,LSFOR
++ LOGICAL LLONG,LSDYN,LSFOR,SD_ROT
++ REAL(q) SD_ROT_Z,SD_ROT_X
+ DIMENSION A(3,3)
+ DIMENSION POSION(3,NIONS)
+ DIMENSION LSFOR (3,NIONS)
+ DIMENSION NITYP(NTYP)
+ CHARACTER*40 FORM
+ CHARACTER*40 SZNAM
++ REAL(q),PARAMETER :: PI=3.1415926536_q
+
+ !-----direct lattice
+ WRITE(IU,'(A40)') SZNAM
+@@ -562,7 +584,14 @@
+ WRITE(IU,FORM) (A(1,I)/SCALE,A(2,I)/SCALE,A(3,I)/SCALE,I=1,3)
+
+ WRITE(IU,'(20I4)') (NITYP(NT),NT=1,NTYP)
+- IF (LSDYN) WRITE(13,'(A18)') 'Selective dynamics'
++ IF (LSDYN) THEN
++ IF (SD_ROT) THEN
++ WRITE(13,'(A30)') 'Transformed selective dynamics'
++ WRITE(13,'(2F)') SD_ROT_Z*180.0/PI,SD_ROT_X*180.0/PI
++ ELSE
++ WRITE(13,'(A18)') 'Selective dynamics'
++ ENDIF
++ ENDIF
+ WRITE(IU,'(A6)')'Direct'
+
+ IF (LSDYN) THEN
+diff -Nur vasp.4.6/poscar.inc vasp.4.6-ct/poscar.inc
+--- vasp.4.6/poscar.inc 2009-07-22 15:05:24.076290013 +0200
++++ vasp.4.6-ct/poscar.inc 2009-09-18 08:01:29.234012116 +0200
+@@ -14,6 +14,9 @@
+ INTEGER NIONS ! actual number of ions
+ INTEGER NIONP ! actual number of ions inc. empty spheres
+ LOGICAL LSDYN ! selective dynamics (yes/ no)
++ LOGICAL SD_ROT ! rotated selective dynamics (yes/ no)
++ REAL(q) SD_ROT_Z ! constraints rotation (z-axis, 1st rot)
++ REAL(q) SD_ROT_X ! constraints rotation (x-axis, 2nd rot)
+ LOGICAL LDIRCO ! positions in direct/recproc. lattice
+ REAL(q), POINTER :: POSION(:,:) ! positions usually same as DYN%POSION
+ LOGICAL,POINTER :: LSFOR(:,:) ! selective dynamics
projects / physik / posic.git / blobdiff