--- /dev/null
+diff -Nur vasp.4.6.orig/dyna.F vasp.4.6-gamma-fullct/dyna.F
+--- vasp.4.6.orig/dyna.F 2009-09-16 10:45:06.050116866 +0200
++++ vasp.4.6-gamma-fullct/dyna.F 2009-11-05 15:09:26.435422082 +0100
+@@ -1164,6 +1164,71 @@
+
+ !*********************************************************************
+ !
++! 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(NI)/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(1)=COS(T_INFO%SD_ROT_Z(NI))*TMP(1)- &
++ & SIN(T_INFO%SD_ROT_Z(NI))*TMP(2)
++ VEC(2)=SIN(T_INFO%SD_ROT_Z(NI))*TMP(1)+ &
++ & COS(T_INFO%SD_ROT_Z(NI))*TMP(2)
++ ENDIF
++ IF (T_INFO%SD_ROT_X(NI)/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(2)=COS(T_INFO%SD_ROT_X(NI))*TMP(2)- &
++ & SIN(T_INFO%SD_ROT_X(NI))*TMP(3)
++ VEC(3)=SIN(T_INFO%SD_ROT_X(NI))*TMP(2)+ &
++ & COS(T_INFO%SD_ROT_X(NI))*TMP(3)
++ ENDIF
++ ELSE
++ IF (T_INFO%SD_ROT_X(NI)/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(2)=COS(-T_INFO%SD_ROT_X(NI))*TMP(2)- &
++ & SIN(-T_INFO%SD_ROT_X(NI))*TMP(3)
++ VEC(3)=SIN(-T_INFO%SD_ROT_X(NI))*TMP(2)+ &
++ & COS(-T_INFO%SD_ROT_X(NI))*TMP(3)
++ ENDIF
++ IF (T_INFO%SD_ROT_Z(NI)/=0._q) THEN
++ DO M=1,3
++ TMP(M)=VEC(M)
++ ENDDO
++ VEC(1)=COS(-T_INFO%SD_ROT_Z(NI))*TMP(1)- &
++ & SIN(-T_INFO%SD_ROT_Z(NI))*TMP(2)
++ VEC(2)=SIN(-T_INFO%SD_ROT_Z(NI))*TMP(1)+ &
++ & COS(-T_INFO%SD_ROT_Z(NI))*TMP(2)
++ ENDIF
++ ENDIF
++ ENDIF
++
++ RETURN
++
++ END
++
++
++!*********************************************************************
++!
+ ! subroutine to set selected force components to zero
+ !
+ !*********************************************************************
+@@ -1184,19 +1249,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 +1312,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.orig/main.F vasp.4.6-gamma-fullct/main.F
+--- vasp.4.6.orig/main.F 2009-09-16 10:45:05.747597375 +0200
++++ vasp.4.6-gamma-fullct/main.F 2009-11-05 15:09:26.449951618 +0100
+@@ -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.orig/pardens.F vasp.4.6-gamma-fullct/pardens.F
+--- vasp.4.6.orig/pardens.F 2009-09-16 10:45:06.294577287 +0200
++++ vasp.4.6-gamma-fullct/pardens.F 2009-11-05 15:09:26.453773187 +0100
+@@ -308,7 +308,7 @@
+ ! I do not call <<write_pard>> here, because I handle the total charge a little bit
+ ! differnt, and the output to <<CHG>> is also differnt
+ do_io CALL OUTPOS(iuchgcar,.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%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ CALL OUTCHG(GRIDC,iuchgcar,.TRUE.,CHTOT)
+
+@@ -320,7 +320,7 @@
+ do_io CLOSE(iuchgcar)
+
+ do_io CALL OUTPOS(iuchg,.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%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+
+ CALL OUTCHG(GRIDC,iuchg,.FALSE.,CHTOT)
+
+@@ -1154,7 +1154,7 @@
+
+ ! and use the usual output routines
+ CALL OUTPOS(iunit,.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%NITYP,.FALSE.,T_INFO%NIONS,DYN%POSION,T_INFO%LSFOR,T_INFO%SD_ROT,T_INFO%SD_ROT_Z,T_INFO%SD_ROT_X)
+ io_end
+
+ CALL OUTCHG(GRIDC,iunit,.FALSE.,CHTOT)
+diff -Nur vasp.4.6.orig/poscar.F vasp.4.6-gamma-fullct/poscar.F
+--- vasp.4.6.orig/poscar.F 2009-09-16 10:45:06.074097686 +0200
++++ vasp.4.6-gamma-fullct/poscar.F 2009-11-05 15:19:56.505775347 +0100
+@@ -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,17 @@
+ 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
++ T_INFO%SD_ROT=.TRUE.
++!-----print some info and convert to rad
++ WRITE(IU6,*) ' rotated selective dynamics switched on'
++!-----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. &
+@@ -288,6 +300,12 @@
+ & DYN%D2C(3,NIOND), &
+ & DYN%VEL(3,NIOND),DYN%D2(3,NIOND),DYN%D3(3,NIOND))
+
++ ALLOCATE(T_INFO%SD_ROT_Z(NIOND))
++ ALLOCATE(T_INFO%SD_ROT_X(NIOND))
++
++ T_INFO%SD_ROT_Z=0._q
++ T_INFO%SD_ROT_X=0._q
++
+ ! alias T_INFO%POSION
+ T_INFO%POSION => DYN%POSION
+
+@@ -298,12 +316,21 @@
+ DYN%D3 =0
+
+ DO NI=1,T_INFO%NIONS
++ IF (T_INFO%SD_ROT) THEN
++ READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI), &
++ & T_INFO%LSFOR(1,NI),T_INFO%LSFOR(2,NI),T_INFO%LSFOR(3,NI), &
++ & T_INFO%SD_ROT_Z(NI), T_INFO%SD_ROT_X(NI)
++!-----convert to rad
++ T_INFO%SD_ROT_Z(NI)=PI*T_INFO%SD_ROT_Z(NI)/180.0
++ T_INFO%SD_ROT_X(NI)=PI*T_INFO%SD_ROT_X(NI)/180.0
++ ELSE
+ IF (T_INFO%LSDYN) THEN
+ READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI), &
+ & T_INFO%LSFOR(1,NI),T_INFO%LSFOR(2,NI),T_INFO%LSFOR(3,NI)
+ ELSE
+ READ(15,*,ERR=400,END=400) DYN%POSION(1,NI),DYN%POSION(2,NI),DYN%POSION(3,NI)
+ ENDIF
++ ENDIF
+ ENDDO
+
+ IF (CSEL=='K') THEN
+@@ -538,17 +565,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), POINTER :: 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,14 +591,26 @@
+ 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'
++ ELSE
++ WRITE(13,'(A18)') 'Selective dynamics'
++ ENDIF
++ ENDIF
+ WRITE(IU,'(A6)')'Direct'
+
+ IF (LSDYN) THEN
+ IF (LLONG) THEN
+ FORM='(3F20.16,3L4)'
++ IF (SD_ROT) THEN
++ FORM='(3F20.16,3L4,2F20.12)'
++ ENDIF
+ ELSE
+ FORM='(3F10.6,3L2)'
++ IF (SD_ROT) THEN
++ FORM='(3F10.6,3L2,2F10.2)'
++ ENDIF
+ ENDIF
+ ELSE
+ IF (LLONG) THEN
+@@ -579,6 +620,12 @@
+ ENDIF
+ ENDIF
+
++ IF(SD_ROT) THEN
++ WRITE(IU,FORM) &
++ & (POSION(1,NI),POSION(2,NI),POSION(3,NI), &
++ & LSFOR(1,NI),LSFOR(2,NI),LSFOR(3,NI), &
++ & SD_ROT_Z(NI)*180.0/PI, SD_ROT_X(NI)*180.0/PI,NI=1,NIONS)
++ ELSE
+ IF (LSDYN) THEN
+ WRITE(IU,FORM) &
+ & (POSION(1,NI),POSION(2,NI),POSION(3,NI), &
+@@ -587,6 +634,7 @@
+ WRITE(IU,FORM) &
+ & (POSION(1,NI),POSION(2,NI),POSION(3,NI),NI=1,NIONS)
+ ENDIF
++ ENDIF
+ IF (.NOT.LLONG) WRITE(IU,*)
+ RETURN
+ END SUBROUTINE
+diff -Nur vasp.4.6.orig/poscar.inc vasp.4.6-gamma-fullct/poscar.inc
+--- vasp.4.6.orig/poscar.inc 2009-09-16 10:45:06.791087379 +0200
++++ vasp.4.6-gamma-fullct/poscar.inc 2009-11-03 11:58:39.388412198 +0100
+@@ -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), POINTER :: SD_ROT_Z(:) ! constraints rotation (z-axis, 1st rot)
++ REAL(q), POINTER :: 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 / commitdiff