C SUBROUTINE ELDEN_ANISO(DEN) IMPLICIT NONE C----------------------------------------------------------------- C--------- SPACE GROUP GENERAL C---This subroutine adds ontribution of atoms with aniso B-value C---to electron desnity. C------------------------------------------------------------------- INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'vitals.fh' INCLUDE 'const.fh' C REAL DEN(N1,N2,N3) REAL GAUSSD(5001) C C---Local variables REAL XA_LIST(3,200),UANISO_LIST(6,200),UANISO_CELL(6), + XYZ_FRAC(3) INTEGER INDSYM_LIST(200) INTEGER I,NX50,NY50,NZ50,IA,IA1,IA11,IL,IL1,NATOM_LIST,ILIST, & IANISO,IG,IDET,ISZL,ISZU,IZ1,IZ,ISYL,ISYU,IY1,IY,IY2, & ISXL,ISXU,IX1,IX,IX2,ID31,IZ2 REAL XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER,SX,SY,SZ,DLIMIT,CSA2,D22, & D238,UEIGEN_MAX,U11,UL,D1,X1,Y1,Z1,XC,YC,ZC,RADZ,SZL,SZU,DZ, & D2,DZO,DZO2,D7,RADY,ZCOSA,SYL,SYU,DZCSA2,DZROX,DZROY,DZCSB, & DY,DYO,DYO2,DYZO,DYZ,DXMIN,DSQMIN,D4,RADX,SXL,SXU,DYZROX, & DXO2PYO2,DX,DXO2,DSQ,DXYO,DXZO,DXYO2,DSQ1, & XDELTA,DXO,DXZO2,DYZO2 REAL CL1(5) REAL FU_ANISO(6,5,MAXNSF),UL_ANISO(6,5),DET_U(5), + XUX(5),XUXZ(5),XUXY(5),UL_INVERT(6,5),FAC,UANISO_F(6), + UANISO_O(6) LOGICAL ERROR CHARACTER LINE*256 C C----Find extension limits for asymmetric unit CALL ASYMLIM_FRAC(DDLIM,XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER) C C---Convert form factors to cell B-values in cell edges. It should be removed C---from here and put to celsym_aniso.fh or something like that CALL FBS2ANISO(NGAUS,MAXNSF,CS_NSFATM,CS_B,FU_ANISO) SX = CS_CELL(1)/NX SY = CS_CELL(2)/NY SZ = CS_CELL(3)/NZ COSAST = (COSA-COSB*COSG)/(SING*SING) DLIMIT = DDLIM C DO I=1,5000 GAUSSD(I)=EXP(-.002*I+.002) ENDDO GAUSSD(5001) = 0.0 CSA2 = COSA*2.0 C --- D22 = SX*SY*SZ*COSZ*SING/TWOPI**1.5 C C----Loop over all atoms NX50 = 500*NX NY50 = 500*NY NZ50 = 500*NZ DO IA=1,N_ATOM IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)-IA1*10 IF(IA11.LE.1) GOTO 500 IF(U_ANISO(2,IA).EQ.0.0.OR.OCCUP(IA).LE.0.0) GOTO 500 C C---If atom has aniso B-value C C---Add NCS here IL = ID_SF(IA) D238 = OCCUP(IA)*D22 IL1 = CS_NELEC(IL) CALL FIND_MAX_EIGEN(U_ANISO(1,IA),UEIGEN_MAX) U11 = 2.0*UEIGEN_MAX IF(IL1.GE.16)U11=AMAX1(U11,0.4) UL = 0.70+U11 IF(IL1.EQ.1)UL = 0.38+U11 D1 = DLIMIT*UL C C---List of all atoms contributing to asymmetric unit. B values of C---symmetry related atoms are different from those of original. CALL MAT2VEC(3,3,CS_ORT_TO_FRAC,XYZ_CRD(1,IA),XYZ_FRAC,ERROR) CALL MAT2VECT(6,6,UUORTH2CELL,U_ANISO(1,IA),UANISO_CELL,ERROR) CALL ATLIST_ANISO(XYZ_FRAC,UANISO_CELL,XA_list,UANISO_LIST, + XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER,INDSYM_LIST,NATOM_LIST) C C----Loop over list of atoms which contribute to "asymmetric" unit DO ILIST=1,NATOM_LIST X1 = XA_list(1,ILIST)*CELL(1) Y1 = XA_list(2,ILIST)*CELL(2) Z1 = XA_list(3,ILIST)*CELL(3) DO IANISO=1,6 UANISO_F(IANISO) = UANISO_LIST(IANISO,ILIST) ENDDO CALL MAT2VECT(6,6,UUCELL2ORTH,UANISO_F,UANISO_O,ERROR) DO IG=1,NGAUS DO IANISO=1,6 UL_ANISO(IANISO,IG) = FU_ANISO(IANISO,IG,IL) + + UANISO_O(IANISO) ENDDO ENDDO C C---Invert B-values and calculate their determinants. CALL ANISO_INVERT(NGAUS,UL_ANISO,DET_U,UL_INVERT) DO IDET = 1,NGAUS IF(DET_U(IDET).LE.0.0) THEN WRITE(LINE,'(A,2I5)')'Problem with atom in ELDEN',IA, + INDSYM_LIST(ILIST) WRITE(*,*)(U_ANISO(IANISO,IA),IANISO=1,6) CALL ERRWRT(-1,LINE) GOTO 500 ENDIF ENDDO DO IG=1,NGAUS CL1(IG) = D238*CS_A(IG,IL)/SQRT(DET_U(IG)) ENDDO XC = X1/SX YC = Y1/SY ZC = Z1/SZ RADZ = SQRT(D1)/(SZ*COSZ) SZL = ZC - RADZ SZU = ZC + RADZ ISZL = INT(SZL+501.0) ISZU = INT(SZU+500.0) DO IZ1 = ISZL,ISZU IZ = IZ1 - 500 IZ2 = MOD(IZ+NZ50,NZ) + 1 IF(IZ2.GT.N3)GO TO 300 IF(IZ2.LT.1) GO TO 300 DZ = IZ*SZ-Z1 D2 = DZ*DZ DZO = DZ*RO_UNIT(3,3) DZO2 = DZO*DZO D7 = D1-DZO2 IF (D7.LT.0.0)GO TO 300 DO IG=1,5 XUXZ(IG) = DZO2*UL_INVERT(3,IG) ENDDO RADY = SQRT(D7)/(SY*SING) ZCOSA = DZ*COSAST/SY SYL = YC - RADY - ZCOSA SYU = YC + RADY - ZCOSA ISYL = INT(SYL+501.0) ISYU = INT(SYU+500.0) DZCSA2 = DZ*CSA2 DZROX = DZ*RO_UNIT(1,3) DZROY = RO_UNIT(2,3)*DZ DZCSB = DZ*COSB DO IY1 = ISYL,ISYU IY = IY1-500 IY2 = MOD(IY+NY50,NY)+1 C---CHECK FOR IY2 WITHIN ASYMMETRIC UNIT IF(IY2.GT.N2)GO TO 290 IF(IY2.LT.1) GO TO 290 DY = IY*SY-Y1 DYO = DY*RO_UNIT(2,2) + DZROY DYO2 = DYO*DYO DYZO = DYO*DZO DYZO2 = 2.0*DYZO DYZ = D2 + (DZCSA2 + DY)*DY DXMIN =-DY*COSG-DZCSB DSQMIN = DYZ-DXMIN**2 DO IG=1,5 XUXY(IG) = XUXZ(IG) + DYO2*UL_INVERT(2,IG) + + DYZO2*UL_INVERT(6,IG) ENDDO D4 = D1-DSQMIN IF(D4.LT.0.0) GOTO 290 RADX = SQRT(D4)/(SX*SING) XDELTA = DXMIN/SX SXL = XC - RADX + XDELTA SXU = XC + RADX + XDELTA ISXL = INT(SXL+501.0) ISXU = INT(SXU+500.0) DYZROX = DY*RO_UNIT(1,2) + DZROX DXO2PYO2 = DYO2 + DZO2 DO IX1 = ISXL,ISXU IX = IX1 - 500 IX2 = MOD(IX+NX50,NX)+1 C---CHECK FOR IX2 WITHIN ASYMMETRIC UNIT IF(IX2.GT.N1)GO TO 280 IF(IX2.LT.1) GO TO 280 DX = IX*SX-X1 DXO = DX*RO_UNIT(1,1) + DYZROX DXO2 = DXO*DXO DSQ = DXO2 + DXO2PYO2 IF(DSQ.GT.D1)GO TO 280 FAC = 0.0 DXYO = DXO*DYO DXZO = DXO*DZO DXYO2 = 2.0*DXYO DXZO2 = 2.0*DXZO DO IG = 1,NGAUS XUX(IG) = XUXY(IG) + UL_INVERT(1,IG)*DXO2 + + UL_INVERT(4,IG)*DXYO2 + + UL_INVERT(5,IG)*DXZO2 DSQ1 = XUX(IG)*250.0+1.0 ID31 = MIN1(DSQ1,5001.1) IF(ID31.LE.5000)FAC = FAC + CL1(IG)*GAUSSD(ID31) ENDDO DEN(IX2,IY2,IZ2) = DEN(IX2,IY2,IZ2)+FAC 280 CONTINUE ENDDO 290 CONTINUE ENDDO 300 CONTINUE ENDDO ENDDO 500 CONTINUE ENDDO RETURN END C SUBROUTINE GRAD_ANISO(DEN,GX,GU1,GQ) C------------------------------------------------------------------ C-------- SPACE GROUP GENEREAL ------------------- C---Calculates gradients. Atoms with anisotropic B-values C------------------------------------------------------------------ INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'vitals.fh' INCLUDE 'const.fh' C REAL GX(*),GU1(*),GQ(*),DEN(N1,N2,N3) C C---Local variables REAL GAUSS(5001) REAL XA_LIST(3,250),FU_ANISO(6,5,MAXNSF),UANISO_INIT(6), + UANISO_LIST(6,250),UL_INVERT(6,5),XUX(5),XUXY(5),XUXZ(5), + UL_ANISO(6,5),DET_U(5),DDET_UDU(6,5), CL1(5),CLB(5),CLQ(5) REAL G1X(5),G1Y(5),G1Z(5),GB0(5),GBXX(5),GBYY(5),GBZZ(5),GBXY(5), + GBXZ(5),GBYZ(5),DUDU(6,6,5),GU(6),GUU(6),XYZ(3),XYZ1(3), + GXP(3),XYZ_FRAC(3),UANISO_F(6),UANISO_O(6),GU11(6) INTEGER INDSYM_LIST(250) LOGICAL ERROR CHARACTER LINE*128 C C-----There are redundant points on the asymmetric unit. Reduce density C----to real asymmetruic unit CALL REDALLI(N1,N2,N3,NX,NY,NZ,ROT,TR,NumSymmetry,DEN) C C----Find extension limits for asymmetric unit CALL ASYMLIM_FRAC(DGLIM,XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER) C C---Convert form factors to cell B-values in cell edges. It should be removed C---from here and put to celsym_aniso.fh or something like that CALL FBS2ANISO(NGAUS,MAXNSF,CS_NSFATM,CS_B,FU_ANISO) SX = CELL(1)/NX SY = CELL(2)/NY SZ = CELL(3)/NZ COSAST = (COSA-COSB*COSG)/(SING*SING) DO I=1,5000 GAUSS(I)=EXP(-.002*I+.002) ENDDO GAUSS(5001)=0.0 C------------------------------------------------------------------- D22=CELL(1)*CELL(2)*CELL(3)*COSZ*SING/TWOPI**1.5 C C---Initialisation should be done before this subroutine IT = 1 IXT = 1 IQ = 1 NX50 = 50*NX NY50 = 50*NY NZ50 = 50*NZ CSA2 = 2.0*COSA DO IA=1,N_ATOM IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)-IA1*10 IF(IA11.LE.0) GOTO 500 IF(IA11.LE.2.OR.U_ANISO(2,IA).EQ.0.0) GOTO 498 C C---Add NCS here D238 = D22*OCCUP(IA) IL = ID_SF(IA) IL1 = CS_NELEC(IL) CALL FIND_MAX_EIGEN(U_ANISO(1,IA),UEIGEN_MAX) U11 = 2.0*UEIGEN_MAX DLIM = DGLIM IF(IL1.EQ.1) DLIM = 3.5 IF(IL1.GE.16)U11 = AMAX1(U11,0.40) RADSQ = DLIM*(U11+0.70) RAD = SQRT(RADSQ) RADZ = RAD/COSZ C C---List of all atoms contributing to asymmetric unit. B values of C---symmetry related atoms are different from those of original. CALL MAT2VEC(3,3,CS_ORT_TO_FRAC,XYZ_CRD(1,IA),XYZ_FRAC,ERROR) CALL MAT2VECT(6,6,UUORTH2CELL,U_ANISO(1,IA),UANISO_INIT,ERROR) CALL ATLIST_ANISO(XYZ_FRAC,UANISO_INIT,XA_list,UANISO_LIST, + XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER,INDSYM_LIST,NATOM_LIST) C C----Loop over list of atoms which contribute to "asymmetric" unit DO ILIST=1,NATOM_LIST C C----Initialisation DO IG=1,NGAUS G1X(IG) = 0.0 G1Y(IG) = 0.0 G1Z(IG) = 0.0 GB0(IG) = 0.0 GBXX(IG) = 0.0 GBYY(IG) = 0.0 GBZZ(IG) = 0.0 GBXY(IG) = 0.0 GBXZ(IG) = 0.0 GBYZ(IG) = 0.0 ENDDO ISYM = INDSYM_LIST(ILIST) X1 = XA_list(1,ILIST)*CELL(1) Y1 = XA_list(2,ILIST)*CELL(2) Z1 = XA_list(3,ILIST)*CELL(3) DO IANISO=1,6 UANISO_F(IANISO) = UANISO_LIST(IANISO,ILIST) ENDDO CALL MAT2VECT(6,6,UUCELL2ORTH,UANISO_F,UANISO_O,ERROR) DO IG=1,NGAUS DO IANISO=1,6 UL_ANISO(IANISO,IG) = FU_ANISO(IANISO,IG,IL) + + UANISO_O(IANISO) ENDDO ENDDO C C---Invert B-values and calculate corresponding derivatives. CALL ANISO_INVERT(NGAUS,UL_ANISO,DET_U,UL_INVERT) DO IDET = 1,NGAUS IF(DET_U(IDET).LE.0.0) THEN WRITE(LINE,'(A,2I5)')'Problem with atom in GRAD ',IA, + INDSYM_LIST(ILIST) CALL ERRWRT(-1,LINE) GOTO 498 ENDIF ENDDO DO IG=1,NGAUS CL1(IG) = D238*CS_A(IG,IL)/SQRT(DET_U(IG)) CLB(IG) = 0.5*CL1(IG) CLQ(IG) = D22*CS_A(IG,IL)/SQRT(DET_U(IG)) ENDDO ZL = Z1 - RADZ ZU = Z1 + RADZ IZL = INT(ZL/SZ+501.) IZU = INT(ZU/SZ+500.) NPOINTS = 0 DO IZ1=IZL,IZU IZ = IZ1-500 IZ2 = MOD(IZ+NZ50,NZ)+1 IF(IZ2.GT.N3)GO TO 300 IF(IZ2.LT.1) GO TO 300 DZ = IZ*SZ-Z1 DZZ = DZ*DZ DZO = DZ*RO_UNIT(3,3) DZO2 = DZO*DZO D7 = RADSQ-DZO2 IF (D7.LT.0.0)GO TO 300 RADY = SQRT(D7)/SING Y1P = Y1 - DZ*COSAST YL = Y1P - RADY YU = Y1P + RADY IYL = INT(YL/SY+501.) IYU = INT(YU/SY+500.) IF(IYU.LT.IYL)GO TO 300 DO IG=1,NGAUS XUXZ(IG) = UL_INVERT(3,IG)*DZO2 ENDDO DZROX = DZ*RO_UNIT(1,3) DZROY = DZ*RO_UNIT(2,3) DZCSA2 = DZ*CSA2 DZCSB = DZ*COSB DO IY1=IYL,IYU IY = IY1-500 IY2 = MOD(IY+NY50,NY)+1 C---CHECK FOR IY2 WITHIN ASYMMETRIC UNIT IF(IY2.GT.N2)GO TO 290 IF(IY2.LT.1) GO TO 290 DY = IY*SY-Y1 DXMIN =-DY*COSG - DZCSB DYZ1 = DZZ + DY*(DZCSA2 + DY) D4 = RADSQ - DYZ1 + DXMIN**2 IF (D4.LT.0.0)GO TO 290 RADX = SQRT(D4)/SING X1P = X1 + DXMIN XL = X1P - RADX XU = X1P + RADX IXL = INT(XL/SX+501.) IXU = INT(XU/SX+500.) IF(IXU.LT.IXL)GO TO 290 DYO = DY*RO_UNIT(2,2) + DZROY DYO2 = DYO*DYO DYZO = DYO*DZO DO IG=1,NGAUS XUXY(IG) = XUXZ(IG) + DYO2*UL_INVERT(2,IG) + + 2.0*UL_INVERT(6,IG)*DYZO ENDDO DYZROX = DY*RO_UNIT(1,2) + DZROX DYO2PZO2 = DYO2 + DZO2 DO IX1=IXL,IXU IX = IX1-500 C---CHECK FOR IX2 WITHIN ASYMMETRIC UNIT IX2 = MOD(IX+NX50,NX)+1 IF(IX2.GT.N1)GO TO 280 IF(IX2.LT.1) GO TO 280 DD = DEN(IX2,IY2,IZ2) IF(DD.EQ.0.0)GOTO 280 DX = IX*SX-X1 DXO = DX*RO_UNIT(1,1) + DYZROX DXO2 = DXO*DXO DSQ = DXO2 + DYO2PZO2 IF(DSQ.GT.RADSQ)GO TO 280 DXYO = DXO*DYO DXZO = DXO*DZO DXYO2 = 2.0*DXYO DXZO2 = 2.0*DXZO C C---We could use here variable integration if we would like to do so. C---But then we could loose speed. (!!!!!) DO IG=1,NGAUS XUX(IG) = XUXY(IG) + DXO2*UL_INVERT(1,IG) + + UL_INVERT(4,IG)*DXYO2 + UL_INVERT(5,IG)*DXZO2 DSQ1 = XUX(IG)*250.0+1.0 ID31 = MIN1(DSQ1,5001.1) IF(ID31.LT.5001) THEN D31 = GAUSS(ID31)*DD G1X (IG) = G1X (IG) + DXO*D31 G1Y (IG) = G1Y (IG) + DYO*D31 G1Z (IG) = G1Z (IG) + DZO*D31 GB0 (IG) = GB0 (IG) + D31 GBXX(IG) = GBXX(IG) + DXO2*D31 GBYY(IG) = GBYY(IG) + DYO2*D31 GBZZ(IG) = GBZZ(IG) + DZO2*D31 GBXY(IG) = GBXY(IG) + DXYO*D31 GBXZ(IG) = GBXZ(IG) + DXZO*D31 GBYZ(IG) = GBYZ(IG) + DYZO*D31 NPOINTS = NPOINTS + 1 ENDIF ENDDO 280 CONTINUE ENDDO 290 CONTINUE ENDDO 300 CONTINUE ENDDO C C---If no grid points was involved in calculation then skip follwoing C---derivative calculations. It might happen for some symmetry related C---atoms or atoms outside asymmetric unit. (In many space groups C---asymmetric unit used here is much larger than real one) IF(NPOINTS.LE.0) GOTO 400 CALL ANISO_INVERT_1DERIV(NGAUS,UL_ANISO,DET_U,UL_INVERT, + DDET_UDU,DUDU) C C----Calculate derivatives wrt current positional parameters GXP(1) = 0.0 GXP(2) = 0.0 GXP(3) = 0.0 DO I=1,6 GU(I) = 0.0 ENDDO GQQ = 0.0 DO IG=1,NGAUS GXP(1) = GXP(1) + CL1(IG)*(UL_INVERT(1,IG)*G1X(IG) + + UL_INVERT(4,IG)*G1Y(IG) + + UL_INVERT(5,IG)*G1Z(IG)) GXP(2) = GXP(2) + CL1(IG)*(UL_INVERT(4,IG)*G1X(IG) + + UL_INVERT(2,IG)*G1Y(IG) + + UL_INVERT(6,IG)*G1Z(IG)) GXP(3) = GXP(3) + CL1(IG)*(UL_INVERT(5,IG)*G1X(IG) + + UL_INVERT(6,IG)*G1Y(IG) + + UL_INVERT(3,IG)*G1Z(IG)) C c---Same thing for B (or) U-values DO I=1,6 GU(I) = GU(I) + CLB(IG)*(DDET_UDU(I,IG)*GB0(IG) + + DUDU(1,I,IG)*GBXX(IG) + DUDU(2,I,IG)*GBYY(IG) + + DUDU(3,I,IG)*GBZZ(IG) + + 2.0*(DUDU(4,I,IG)*GBXY(IG) + + DUDU(5,I,IG)*GBXZ(IG) + DUDU(6,I,IG)*GBYZ(IG))) ENDDO GQQ = GQQ + CLQ(IG)*GB0(IG) ENDDO C C---If current atom is symmetry related to original atom then calculate C---derivative of original one by using symmetry operators (simple chain C---rule) IF(ISYM.GT.1) THEN CALL MAT2VECT(3,3,CS_FRAC_TO_ORT,GXP,XYZ,ERROR) CALL MAT2VECT(3,3,ROT(1,1,ISYM),XYZ,XYZ1,ERROR) CALL MAT2VECT(3,3,CS_ORT_TO_FRAC,XYZ1,GXP,ERROR) CALL MAT2VEC(6,6,UUCELL2ORTH,GU,GU11,ERROR) CALL MAT2VEC(6,6,RealSymm_Aniso(1,1,ISYM),GU11,GUU,ERROR) CALL MAT2VEC(6,6,UUORTH2CELL,GUU,GU,ERROR) ENDIF GX(IXT ) = GX(IXT ) + GXP(1) GX(IXT+1) = GX(IXT+1) + GXP(2) GX(IXT+2) = GX(IXT+2) + GXP(3) DO IANISO = 0,5 GU1(IT+IANISO) = GU1(IT+IANISO) - GU(IANISO+1) ENDDO cd IF(OCCUP_REF(IA).GT.1) GQ(IQ) = GQ(IQ) + GQQ GQ(IQ) = GQ(IQ) + GQQ 400 CONTINUE ENDDO C 498 CONTINUE IXT = IXT + 3 IF(U_ANISO(2,IA).EQ.0.0) THEN IT = IT + 1 ELSE IT = IT + 6 ENDIF IQ = IQ + 1 500 CONTINUE ENDDO RETURN END C SUBROUTINE ANISO_INVERT(NGAUS,BL_ANISO,DET_B,BL_INVERT) C C---Calculates inversion of matrices and its determinant (of original matrix) C---for anisotropic C---B-values. Uses symmetricity of Aniso B tensor. This subrotine uses C---Kronekker's rule ofr calculation inversion of matrices REAL BL_ANISO(6,NGAUS),DET_B(NGAUS),BL_INVERT(6,NGAUS) C C---First calculate minors of anisotropic tensor DO IG=1,NGAUS BL_INVERT(1,IG) = + BL_ANISO(2,IG)*BL_ANISO(3,IG) - BL_ANISO(6,IG)**2 BL_INVERT(2,IG) = + BL_ANISO(1,IG)*BL_ANISO(3,IG) - BL_ANISO(5,IG)**2 BL_INVERT(3,IG) = + BL_ANISO(1,IG)*BL_ANISO(2,IG) - BL_ANISO(4,IG)**2 BL_INVERT(4,IG) = + -BL_ANISO(3,IG)*BL_ANISO(4,IG) + + BL_ANISO(5,IG)*BL_ANISO(6,IG) BL_INVERT(5,IG) = + BL_ANISO(4,IG)*BL_ANISO(6,IG) - + BL_ANISO(5,IG)*BL_ANISO(2,IG) BL_INVERT(6,IG) = + -BL_ANISO(1,IG)*BL_ANISO(6,IG) + + BL_ANISO(4,IG)*BL_ANISO(5,IG) DET_B(IG) = 0.0 ENDDO C C----Calulate determinants. DO IG=1,NGAUS DET_B(IG) = (BL_ANISO(1,IG)*BL_INVERT(1,IG)+ + BL_ANISO(4,IG)*BL_INVERT(4,IG)+ + BL_ANISO(5,IG)*BL_INVERT(5,IG)) ENDDO C C---Now calculate inversion of matrices using Kronekkers rule DO IG=1,NGAUS DO IANISO = 1,6 BL_INVERT(IANISO,IG) = BL_INVERT(IANISO,IG)/DET_B(IG) ENDDO ENDDO C C----Derivatives here or not here that is the question END C SUBROUTINE ANISO_INVERT_1DERIV(NGAUS,UL_ANISO,DET_U,UL_INVERT, + DDET_UDU,DUDU) C C--Calculates derivatives of inverted matrix with respect to matrix elements C--necessary to calculate gradients REAL UL_ANISO(6,NGAUS),DET_U(NGAUS),UL_INVERT(6,NGAUS) REAL DDET_UDU(6,NGAUS),DUDU(6,6,NGAUS) DO IG=1,NGAUS DDET_UDU(1,IG) = UL_INVERT(1,IG) DDET_UDU(2,IG) = UL_INVERT(2,IG) DDET_UDU(3,IG) = UL_INVERT(3,IG) DDET_UDU(4,IG) = 2.0*UL_INVERT(4,IG) DDET_UDU(5,IG) = 2.0*UL_INVERT(5,IG) DDET_UDU(6,IG) = 2.0*UL_INVERT(6,IG) DO I=1,6 DO J=1,6 DUDU(I,J,IG) = -UL_INVERT(I,IG)*DDET_UDU(J,IG) ENDDO ENDDO DUDU(2,1,IG) = DUDU(2,1,IG) + UL_ANISO(3,IG)/DET_U(IG) DUDU(3,1,IG) = DUDU(3,1,IG) + UL_ANISO(2,IG)/DET_U(IG) DUDU(6,1,IG) = DUDU(6,1,IG) - UL_ANISO(6,IG)/DET_U(IG) DUDU(1,2,IG) = DUDU(1,2,IG) + UL_ANISO(3,IG)/DET_U(IG) DUDU(3,2,IG) = DUDU(3,2,IG) + UL_ANISO(1,IG)/DET_U(IG) DUDU(5,2,IG) = DUDU(5,2,IG) - UL_ANISO(5,IG)/DET_U(IG) DUDU(1,3,IG) = DUDU(1,3,IG) + UL_ANISO(2,IG)/DET_U(IG) DUDU(2,3,IG) = DUDU(2,3,IG) + UL_ANISO(1,IG)/DET_U(IG) DUDU(4,3,IG) = DUDU(4,3,IG) - UL_ANISO(4,IG)/DET_U(IG) DUDU(3,4,IG) = DUDU(3,4,IG) - 2.0*UL_ANISO(4,IG)/DET_U(IG) DUDU(4,4,IG) = DUDU(4,4,IG) - UL_ANISO(3,IG)/DET_U(IG) DUDU(5,4,IG) = DUDU(5,4,IG) + UL_ANISO(6,IG)/DET_U(IG) DUDU(6,4,IG) = DUDU(6,4,IG) + UL_ANISO(5,IG)/DET_U(IG) DUDU(2,5,IG) = DUDU(2,5,IG) - 2.0*UL_ANISO(5,IG)/DET_U(IG) DUDU(4,5,IG) = DUDU(4,5,IG) + UL_ANISO(6,IG)/DET_U(IG) DUDU(5,5,IG) = DUDU(5,5,IG) - UL_ANISO(2,IG)/DET_U(IG) DUDU(6,5,IG) = DUDU(6,5,IG) + UL_ANISO(4,IG)/DET_U(IG) DUDU(1,6,IG) = DUDU(1,6,IG) - 2.0*UL_ANISO(6,IG)/DET_U(IG) DUDU(4,6,IG) = DUDU(4,6,IG) + UL_ANISO(5,IG)/DET_U(IG) DUDU(5,6,IG) = DUDU(5,6,IG) + UL_ANISO(4,IG)/DET_U(IG) DUDU(6,6,IG) = DUDU(6,6,IG) - UL_ANISO(1,IG)/DET_U(IG) ENDDO END C SUBROUTINE FBS2ANISO(NGAUS,MXSCAT,NSCAT,FBB,FU_ANISO) C IMPLICIT NONE INCLUDE 'const.fh' C C---Converts B-values of atomic form factors to anisotropic form. C---It is necessary if cell is not orthogonal INTEGER NGAUS,MXSCAT,NSCAT REAL FBB(NGAUS,MXSCAT) REAL FU_ANISO(6,NGAUS,MXSCAT) C INTEGER ISCAT,IG DO ISCAT=1,NSCAT DO IG=1,NGAUS FU_ANISO(1,IG,ISCAT) = FBB(IG,ISCAT)/PISQ8 FU_ANISO(2,IG,ISCAT) = FBB(IG,ISCAT)/PISQ8 FU_ANISO(3,IG,ISCAT) = FBB(IG,ISCAT)/PISQ8 FU_ANISO(4,IG,ISCAT) = 0.0 FU_ANISO(5,IG,ISCAT) = 0.0 FU_ANISO(6,IG,ISCAT) = 0.0 ENDDO ENDDO RETURN END C SUBROUTINE BISO2ANISO(BISO,UANISO) INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'const.fh' C C---Converts B-values of atomic form factors to anisotropic U form. REAL UANISO(6) UANISO(1) = BISO UANISO(2) = BISO UANISO(3) = BISO UANISO(4) = 0.0 UANISO(5) = 0.0 UANISO(6) = 0.0 RETURN END C SUBROUTINE UANISO_CELL2ORTH(U_CELL_ANISO,U_ORTH_ANISO) C C----Converts anisotropic U values from cell edges to orthogonal system C----Im REAL INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' REAL U_CELL_ANISO(6),U_ORTH_ANISO(6) DO I=1,6 U_ORTH_ANISO(I) = 0.0 DO J=1,6 U_ORTH_ANISO(I) = U_ORTH_ANISO(I)+ + U_CELL_ANISO(J) * UUCELL2ORTH(J,I) ENDDO ENDDO RETURN END SUBROUTINE UANISO_ORTH2CELL(U_CELL_ANISO,U_ORTH_ANISO) C C----Converts anisotropic U values from orhtogonal system to cell edges C----In REAL INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' REAL U_CELL_ANISO(6),U_ORTH_ANISO(6) DO I=1,6 U_CELL_ANISO(I) = 0.0 DO J=1,6 U_CELL_ANISO(I) = U_CELL_ANISO(I) + + U_ORTH_ANISO(J)*UUORTH2CELL(J,I) ENDDO ENDDO RETURN END C SUBROUTINE ATLIST_ANISO(XFRAC_input, + UANISO_INIT,XA_list,UANISO_LIST, + XLOW,YLOW,ZLOW,XUPPER,YUPPER,ZUPPER,INDSYM_list, + NATOM_list) INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'vitals.fh' C C----Finds list of atoms which has something to do with enlarged C----asymmetric unit after application of symmetry operations. C----Could be used for atoms with anisotropic B-values. REAL XFRAC_input(3),XA_list(3,*),UANISO_INIT(6), + UANISO_LIST(6,*) INTEGER INDSYM_list(*) C---------------------------------------------------------------- C-----Input C----- C----- XFRAC_input, C----- YFRAC_input, C----- ZFRAC_input positional paramters of atoms in fractional C----- coordinates C----- UANISO_INIT thermal parameters in cell edges (angstroms^2) C----- XLOW, C----- YLOW, C----- ZLOW, C----- XUPPER, C----- YUPPER, C----- ZUPPER lower and upper limits of asymmetric unit for C----- enlarging (in fractional coordiates) C----- C-----Output C----- C----- Xa_list list of coordinates C----- UANISO_LIST U-values for symmetry related atoms C----- INDSYM_list corresponding symmetry operation C----- NATOM_list number of atoms in the list C----- C--------------------------------------------------------------- NATOM_list = 0 DO ISYM=1,NumSymmetry XT = RealSymmMatrx(1,1,ISYM)*XFRAC_input(1) + + RealSymmMatrx(1,2,ISYM)*XFRAC_input(2) + + RealSymmMatrx(1,3,ISYM)*XFRAC_input(3) + + RealSymmMatrx(1,4,ISYM) YT = RealSymmMatrx(2,1,ISYM)*XFRAC_input(1) + + RealSymmMatrx(2,2,ISYM)*XFRAC_input(2) + + RealSymmMatrx(2,3,ISYM)*XFRAC_input(3) + + RealSymmMatrx(2,4,ISYM) ZT = RealSymmMatrx(3,1,ISYM)*XFRAC_input(1) + + RealSymmMatrx(3,2,ISYM)*XFRAC_input(2) + + RealSymmMatrx(3,3,ISYM)*XFRAC_input(3) + + RealSymmMatrx(3,4,ISYM) C C---Check if atom is inside enlarged asymmetric unit 20 IF(XT.LT.XLOW)GO TO 21 XT=XT-1. GO TO 20 21 XT=XT+1. IF(XT.LT.XLOW) GO TO 21 30 IF(YT.LT.YLOW) GO TO 31 YT=YT-1. GO TO 30 31 YT=YT+1. IF(YT.LT.YLOW)GO TO 31 40 IF(ZT.LT.ZLOW)GO TO 41 ZT=ZT-1. GO TO 40 41 ZT=ZT+1. IF(ZT.LT.ZLOW)GO TO 41 70 IF(XT.GT.XUPPER)GO TO 90 IF(YT.GT.YUPPER)GO TO 90 IF(ZT.GT.ZUPPER)GO TO 90 C C---Atom is inside enlarged asymmetric unit. Store its coordinates and C---symmetry operation by which this atom came into asymmetric unit NATOM_list = NATOM_list + 1 XA_list(1,NATOM_list) = XT XA_list(2,NATOM_list) = YT XA_list(3,NATOM_list) = ZT INDSYM_list(NATOM_list) = ISYM DO I=1,6 UANISO_LIST(I,NATOM_list) = 0.0 DO J=1,6 UANISO_LIST(I,NATOM_list) = UANISO_LIST(I,NATOM_list) + + RealSymm_Aniso(J,I,ISYM)*UANISO_INIT(J) ENDDO ENDDO 90 CONTINUE ENDDO RETURN END C SUBROUTINE FIND_CONV_MATRIX(R_IN,R_OUT) C C----Calculates conversion matrix for anisotropic u corresbonding C----to coordinate conversion matrix R_in REAL R_IN(3,3),R_OUT(6,6) C INTEGER II,JJ DO II=1,3 DO JJ = 1,3 R_OUT(II,JJ) = R_IN(JJ,II)*R_IN(JJ,II) ENDDO R_OUT(4,II) = 2.0*R_IN(II,1)*R_IN(II,2) R_OUT(5,II) = 2.0*R_IN(II,1)*R_IN(II,3) R_OUT(6,II) = 2.0*R_IN(II,2)*R_IN(II,3) R_OUT(II,4) = R_IN(1,II)*R_IN(2,II) R_OUT(II,5) = R_IN(1,II)*R_IN(3,II) R_OUT(II,6) = R_IN(2,II)*R_IN(3,II) ENDDO R_OUT(4,4) = R_IN(1,1)*R_IN(2,2) + R_IN(1,2)*R_IN(2,1) R_OUT(4,5) = R_IN(1,1)*R_IN(3,2) + R_IN(1,2)*R_IN(3,1) R_OUT(4,6) = R_IN(2,1)*R_IN(3,2) + R_IN(2,2)*R_IN(3,1) R_OUT(5,4) = R_IN(1,1)*R_IN(2,3) + R_IN(1,3)*R_IN(2,1) R_OUT(5,5) = R_IN(1,1)*R_IN(3,3) + R_IN(1,3)*R_IN(3,1) R_OUT(5,6) = R_IN(2,1)*R_IN(3,3) + R_IN(2,3)*R_IN(3,1) R_OUT(6,4) = R_IN(1,2)*R_IN(2,3) + R_IN(1,3)*R_IN(2,2) R_OUT(6,5) = R_IN(1,2)*R_IN(3,3) + R_IN(1,3)*R_IN(3,2) R_OUT(6,6) = R_IN(2,2)*R_IN(3,3) + R_IN(2,3)*R_IN(3,2) C RETURN END C SUBROUTINE APPLY_CONV2U(UANISO,UANISO1,R_CONV) C C---Applies conversion matrix to aniso U REAL UANISO(6),UANISO1(6),R_CONV(6,6) LOGICAL ERROR C CALL MAT2VECT(6,6,R_CONV,UANISO,UANISO1,ERROR) RETURN END C SUBROUTINE ALLBS2ANISOU IMPLICIT NONE C C---Converts isotropic Bs to anisotropic U values INCLUDE 'vitals.fh' INCLUDE 'atom_com.fh' INTEGER IA C C----If atom is not already anisotropic make it so. NANISO_ATOMS = 0 N_ANISO = 0 DO IA=1,N_ATOM IF(U_ANISO(2,IA).LE.0.0) THEN IF(CS_ELEMENT(ID_SF(IA)).NE.'H'.AND. & CS_ELEMENT(ID_SF(IA)).NE.'H-1') THEN U_ANISO(2,IA) = U_ANISO(1,IA) U_ANISO(3,IA) = U_ANISO(1,IA) U_ANISO(4,IA) = 0.0 U_ANISO(5,IA) = 0.0 U_ANISO(6,IA) = 0.0 NANISO_ATOMS = NANISO_ATOMS + 1 N_ANISO = N_ANISO + 1 ENDIF ENDIF ENDDO cd NANISO_ATOMS = N_ATOM cd N_ANISO = N_ATOM END C SUBROUTINE CHECK_TLS_TENSORTS(TMAT,LMAT,SMAT,IFLAG) IMPLICIT NONE C INCLUDE 'const.fh' INTEGER IFLAG REAL TMAT(6),LMAT(6),SMAT(8) REAL UANISO0(6,6),EVALUE(6),WORKSPACE(100) INTEGER I,J,LWORK,INFO C IFLAG = 0 LWORK = 10 UANISO0(1,1) = TMAT(1) UANISO0(2,2) = TMAT(2) UANISO0(3,3) = TMAT(3) UANISO0(1,2) = TMAT(4) UANISO0(1,3) = TMAT(5) UANISO0(2,3) = TMAT(6) UANISO0(4,4) = LMAT(1) UANISO0(5,5) = LMAT(2) UANISO0(6,6) = LMAT(3) UANISO0(4,5) = LMAT(4) UANISO0(4,6) = LMAT(5) UANISO0(5,6) = LMAT(6) C C---Now add S tensor. Assume s11 + s22 + s33 = 0.0 UANISO0(1,4) = (SMAT(2)-SMAT(1))/3.0 UANISO0(2,5) = (SMAT(2) + 2.0*SMAT(1))/3.0 UANISO0(3,5) = -(SMAT(1) + 2.0*SMAT(2))/3.0 UANISO0(2,4) = SMAT(3) UANISO0(3,4) = SMAT(4) UANISO0(3,5) = SMAT(5) UANISO0(1,5) = SMAT(6) UANISO0(2,6) = SMAT(7) UANISO0(3,6) = SMAT(8) DO I=1,5 DO J=I+1,6 UANISO0(J,I) = UANISO0(I,J) ENDDO ENDDO CALL SSYEV('V','U',3,UANISO0,3,EVALUE,WORKSPACE,LWORK,INFO) DO I=1,6 IF(EVALUE(I).LT.0.0) THEN WRITE(*,*)'Eigenvalue is less than 0. Reset it' EVALUE(I) = 0.0 IFLAG = 1 ENDIF IF(EVALUE(I).GT.20.0/(RTODEG*RTODEG)) THEN WRITE(*,*)'Eigenvalue too large. Reset it' EVALUE(I) = 20.0/(RTODEG*RTODEG) IFLAG = 1 ENDIF cd EVALUE(I) = AMAX1(0.0,AMIN1(EVALUE(I),30.0/(RTODEG*RTODEG))) ENDDO cd CALL EIGEN2U(EVALUE,UANISO0,LMAT) C C---Now return back to TLS parameters RETURN END C SUBROUTINE CHECK_U_VALUES IMPLICIT NONE C C---Checks validity of U values. In isotropic case if U is C---outside predefined limit UMIN, UMAX makes it inside it. In case C---anisotropic U finds eigenvalues checks if they are outside C---predefined limit. If they are makes them inside INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'vitals.fh' INCLUDE 'refi_flags.fh' INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' C REAL UANISO0(3,3),UANISO1(6),EVALUE(3),WORKSPACE(20) CHARACTER LINE*256,AT_FULL*15 C INTEGER IA,LWORK,INFO,I,II,IERROR,IANISO LWORK = 20 DO IA=1,N_ATOM IF(U_ANISO(2,IA).EQ.0.0) THEN C C---Atom is isotropic. put it inside predefined limit. IF(U_ANISO(1,IA).LT.BResetMin) THEN U_ANISO(1,IA) = BResetMin ENDIF IF(U_ANISO(1,IA).GT.BResetMax) THEN U_ANISO(1,IA) = BResetMax ENDIF ELSE C C---Atom is anisotropic UANISO0(1,1) = U_ANISO(1,IA) UANISO0(2,2) = U_ANISO(2,IA) UANISO0(3,3) = U_ANISO(3,IA) UANISO0(1,2) = U_ANISO(4,IA) UANISO0(1,3) = U_ANISO(5,IA) UANISO0(2,3) = U_ANISO(6,IA) UANISO0(2,1) = UANISO0(1,2) UANISO0(3,1) = UANISO0(1,3) UANISO0(3,2) = UANISO0(2,3) CALL SSYEV('V','U',3,UANISO0,3,EVALUE,WORKSPACE,LWORK,INFO) C C---Check if atom has small eigen value IF(EVALUE(1).LT.BResetMin) THEN CALL FULL_ATOM_NAME(IA,AT_FULL,IERROR) WRITE(LINE,'(A,A,A)')'Atom ',AT_FULL, + ' tends to have small eigenvalue. Reset it' CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,3F10.4)')'Eigenvalues ', + EVALUE(1),EVALUE(2),EVALUE(3) CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,6F10.4)')' U before setting ', + (U_ANISO(II,IA),II=1,6) CALL ERRWRT(-1,LINE) EVALUE(1) = AMAX1(BResetMin,EVALUE(1)) EVALUE(2) = AMAX1(BResetMin,EVALUE(2)) EVALUE(3) = AMAX1(BResetMin,EVALUE(3)) CALL EIGEN2U(EVALUE,UANISO0,UANISO1) DO IANISO=1,6 U_ANISO(IANISO,IA) = UANISO1(IANISO) ENDDO WRITE(LINE,'(A,6F10.4)')' U after setting ', + (U_ANISO(II,IA),II=1,6) CALL ERRWRT(-1,LINE) ENDIF C C---Check if atom has too big eigenvalue IF(EVALUE(3).GT.BResetMax ) THEN CALL FULL_ATOM_NAME(IA,AT_FULL,IERROR) WRITE(LINE,'(A,A,A)')'Atom ',AT_FULL, + ' tends to have large eigenvalue. Reset it' CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,3F10.3)')'Eigenvalues ', + EVALUE(1),EVALUE(2),EVALUE(3) CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,6F10.4)')' U before setting ', + (U_ANISO(II,IA),II=1,6) EVALUE(1) = AMIN1(BResetMax,EVALUE(1)) EVALUE(2) = AMIN1(BResetMax,EVALUE(2)) EVALUE(3) = AMIN1(BResetMax,EVALUE(3)) CALL EIGEN2U(EVALUE,UANISO0,UANISO1) DO IANISO=1,6 U_ANISO(IANISO,IA) = UANISO1(IANISO) ENDDO WRITE(LINE,'(A,6F10.4)')' U after setting ', + (U_ANISO(II,IA),II=1,6) CALL ERRWRT(-1,LINE) ENDIF ENDIF ENDDO END C SUBROUTINE EIGEN2U(U0,EV,U1) IMPLICIT NONE C C----Converts eigenvalues and eigenvaectors of U tensor to U tensor C----itself. eigenvalues are U0(1),U0(2),U0(2), Evector is eigenvector C----first eigenvalue and so on. C----U1 is aniso U tensor. First 3 elements are diagonal terms C----recall that U = EVt U0diag EV REAL U0(3),U1(6),EV(3,3) C C---Local variables REAL EV1(3,3),EV2(3,3) INTEGER I,J C LOGICAL ERROR C DO I=1,3 DO J=1,3 EV1(I,J) = EV(I,J)*U0(J) ENDDO ENDDO CALL MAT2MATT(3,3,EV1,EV,EV2,ERROR) U1(1) = EV2(1,1) U1(2) = EV2(2,2) U1(3) = EV2(3,3) U1(4) = EV2(1,2) U1(5) = EV2(1,3) U1(6) = EV2(2,3) RETURN END C SUBROUTINE FIND_MIN_EIGEN(UANISO_INIT,UEIGEN_MIN) IMPLICIT NONE C C----Find minimum eigen value for aniso u value. Aniso U values C----are assumed to be in orthogonal coordinates REAL UANISO_INIT(6) REAL UEIGEN_MIN C REAL UANISO0(3,3),EVALUE(3),WORKSPACE(20) INTEGER LWORK,INFO C C--First copy u values for special storage matrix. LWORK = 20 UANISO0(1,1) = UANISO_INIT(1) UANISO0(2,2) = UANISO_INIT(2) UANISO0(3,3) = UANISO_INIT(3) UANISO0(1,2) = UANISO_INIT(4) UANISO0(1,3) = UANISO_INIT(5) UANISO0(2,3) = UANISO_INIT(6) UANISO0(2,1) = UANISO0(1,2) UANISO0(3,1) = UANISO0(1,3) UANISO0(3,2) = UANISO0(2,3) CALL SSYEV('N','U',3,UANISO0,3,EVALUE,WORKSPACE,LWORK,INFO) C C---Find maximum of eigenvalues UEIGEN_MIN = EVALUE(1) END C SUBROUTINE FIND_MAX_EIGEN(UANISO_INIT,UEIGEN_MAX) IMPLICIT NONE C C----Find maximum eigen value for aniso u value. Aniso U values C----are assumed to be in orthogonal coordinates REAL UANISO_INIT(6) REAL UEIGEN_MAX C REAL UANISO0(3,3),EVALUE(3),WORKSPACE(20) INTEGER LWORK,INFO C C--First copy u values for special storage matrix. LWORK = 20 UANISO0(1,1) = UANISO_INIT(1) UANISO0(2,2) = UANISO_INIT(2) UANISO0(3,3) = UANISO_INIT(3) UANISO0(1,2) = UANISO_INIT(4) UANISO0(1,3) = UANISO_INIT(5) UANISO0(2,3) = UANISO_INIT(6) UANISO0(2,1) = UANISO0(1,2) UANISO0(3,1) = UANISO0(1,3) UANISO0(3,2) = UANISO0(2,3) CALL SSYEV('N','U',3,UANISO0,3,EVALUE,WORKSPACE,LWORK,INFO) C C---Find maximum of eigenvalues UEIGEN_MAX = EVALUE(3) END