From 839a6628d140c5c32c108d63ce59b4c7fa2b8858 Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 8 Feb 2010 14:24:12 +0100 Subject: [PATCH] added the all stoms crt patch --- vasp_tools/sd_rot_all-atoms.patch | 407 ++++++++++++++++++++++++++++++ 1 file changed, 407 insertions(+) create mode 100644 vasp_tools/sd_rot_all-atoms.patch diff --git a/vasp_tools/sd_rot_all-atoms.patch b/vasp_tools/sd_rot_all-atoms.patch new file mode 100644 index 0000000..548f3c6 --- /dev/null +++ b/vasp_tools/sd_rot_all-atoms.patch @@ -0,0 +1,407 @@ +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 <> here, because I handle the total charge a little bit + ! differnt, and the output to <> 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 -- 2.39.2