From 197db0fb9df39ac19c8c09969476a8765fad2bc7 Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 21 Sep 2009 19:32:04 +0200 Subject: [PATCH] patch to enable rotation of the constraints basis (test in progress!) --- vasp_tools/sd_rot.patch | 407 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 407 insertions(+) create mode 100644 vasp_tools/sd_rot.patch diff --git a/vasp_tools/sd_rot.patch b/vasp_tools/sd_rot.patch new file mode 100644 index 0000000..c0c15f2 --- /dev/null +++ b/vasp_tools/sd_rot.patch @@ -0,0 +1,407 @@ +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 -- 2.20.1