C C $Id: getxyz.F,v 1.2 1998/07/16 16:39:58 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE GETXYZ1(NATOMS1,ZMAT1,IZMAT1,IFIRST,XYZ1,IERROR) C C COMPUTES CARTESIAN COORDINATES FROM INTERNAL COORDINATES. C C INPUT: C ------ C NATOMS1 = NUMBER OF ATOMS. C ZMAT1 = BOND, ANGLE, DIHEDRAL ANGLE VALUES. UNITS SHOULD BE C ANGSTROMS, RADIANS, RADIANS, RESPECTIVELY. C IZMAT1 = REFERENCE ATOMS DEFINING ZMAT1. C IFIRST = FIRST ATOM FOR WHICH CARTESIAN COORDINATES ARE NEEDED. C THIS IS NORMALLY EQUAL TO 1, BUT WHEN COMPUTING JACOBIAN C MATRIX IN GRADIENT CALCULATION, THIS NUMBER CAN BE C HIGHER SO THAT TIME IS NOT WASTED COMPUTING COORDINATES C THAT DO NOT CHANGE. C C RETURNED: C --------- C XYZ1 = CARTESIAN COORDINATES. C IERROR = ERROR CODE: IERROR=0 --> SUCCESSFUL CALCULATION. C IERROR=1 --> THREE ATOMS IN A STRAIGHT LINE C PREVENTED ACCURATE DETERMINATION C OF CARTESIAN COORDINATES. C C -- S. DIXON C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" DIMENSION ZMAT1(3,*),IZMAT1(3,*),XYZ1(3,*) DIMENSION P(3,3),V(3) IERROR = 0 C C COORDINATES FOR ATOMS 1, 2, AND 3 ARE HANDLED SEPARATELY. C XYZ1(1,1) = 0.0D0 XYZ1(2,1) = 0.0D0 XYZ1(3,1) = 0.0D0 IF(NATOMS1.EQ.1) GO TO 1000 BLENG = ZMAT1(1,2) XYZ1(1,2) = BLENG XYZ1(2,2) = 0.0D0 XYZ1(3,2) = 0.0D0 IF(NATOMS1.EQ.2) GO TO 1000 BLENG = ZMAT1(1,3) ANGLE = ZMAT1(2,3) NA = IZMAT1(1,3) NB = IZMAT1(2,3) A3 = ANGLE - 90.0D0*RADIAN XREF = XYZ1(1,2) C C THE 3-1-2 NUMBERING SCHEME HAS A DIFFERENT CONVENTION: C IF(NA.EQ.1.AND.NB.EQ.2)THEN A3 = -A3 XREF = XYZ1(1,1) ENDIF C XYZ1(1,3) = XREF + BLENG*SIN(A3) XYZ1(2,3) = BLENG*COS(A3) XYZ1(3,3) = 0.0D0 IF(NATOMS1.EQ.3) GO TO 1000 C C NOW BEGIN GENERAL PROCEDURE FOR DETERMINING CARTESIAN COORDINATES C FROM INTERNAL COORDINATES. C IF(IFIRST.GT.4)THEN ISTART = IFIRST ELSE ISTART = 4 ENDIF DO 500 I=ISTART,NATOMS1 NA = IZMAT1(1,I) NB = IZMAT1(2,I) BLENG = ZMAT1(1,I) ANGLE = ZMAT1(2,I) XBA = XYZ1(1,NA) - XYZ1(1,NB) YBA = XYZ1(2,NA) - XYZ1(2,NB) ZBA = XYZ1(3,NA) - XYZ1(3,NB) RBASQR = XBA*XBA + YBA*YBA + ZBA*ZBA RBA = DSQRT(RBASQR) C C IF THE I-NA-NB BOND ANGLE IS CLOSE TO 0 OR 180 DEGREES, THEN JUST C PUT ATOM I IN LINE WITH ATOMS NA AND NB. THERE COULD BE PROBLEMS C LATER, HOWEVER. C ANGDEG = ANGLE*DEGREE DIFF = ABS(ABS(ANGDEG) - 180.0D0) IF(DIFF.LT.0.1D0)THEN SCALE = BLENG/RBA IF(ABS(ANGDEG).LT.1.0D0) SCALE = -SCALE XYZ1(1,I) = XYZ1(1,NA) + SCALE*XBA XYZ1(2,I) = XYZ1(2,NA) + SCALE*YBA XYZ1(3,I) = XYZ1(3,NA) + SCALE*ZBA GO TO 500 ENDIF NC = IZMAT1(3,I) XBC = XYZ1(1,NC) - XYZ1(1,NB) YBC = XYZ1(2,NC) - XYZ1(2,NB) ZBC = XYZ1(3,NC) - XYZ1(3,NB) RBCSQR = XBC*XBC + YBC*YBC + ZBC*ZBC C C AT THIS POINT WE KNOW ATOMS I, NA, AND NB DO NOT LIE IN A C STRAIGHT LINE. HOWEVER, IF THE NA-NB-NC ANGLE IS 0 OR 180 C DEGREES, THEN THE COORDINATES OF ATOM I CANNOT BE COMPUTED C WITH CONFIDENCE. C CABC = (XBA*XBC + YBA*YBC + ZBA*ZBC)/(RBASQR*RBCSQR) CDIFF = ABS(ABS(CABC) - 1.0D0) IF(CDIFF.LT.1.0D-6)THEN IERROR = 1 IF(MYID.EQ.0) THEN C-RDC WRITE(IUNIT(2),20) I,NA,NB,NC ELSE C-RDC WRITE(IUNIT(12),20) I,NA,NB,NC ENDIF 20 FORMAT(/' COORDINATES FOR ATOM ',I5,' CANNOT BE COMPUTED', . ' WITH CONFIDENCE BECAUSE'/' ATOMS ',I5,', ',I5,',', . ' AND ',I5,' LIE IN A STRAIGHT LINE') GO TO 1000 ENDIF C C CONSTRUCT A HOUSEHOLDER TRANSFORMATION P THAT WOULD BRING THE C VECTOR (XBA,YBA,ZBA) IN LINE WITH THE POSITIVE X AXIS. NO C TRANSFORMATION IS NECESSARY IF THIS VECTOR IS A POSITIVE C MULTIPLE OF (1,0,0). C YZBA = YBA**2 + ZBA**2 IF(XBA.LT.0.0D0.OR.YZBA.GT.1.0D-12)THEN ITRAN = 1 V(1) = XBA - RBA V(2) = YBA V(3) = ZBA BETAH = 2.0D0/(V(1)**2 + V(2)**2 + V(3)**2) DO 40 J=1,3 DO 30 K=1,3 P(J,K) = -BETAH*V(J)*V(K) 30 CONTINUE P(J,J) = P(J,J) + 1.0D0 40 CONTINUE C C TRANSFORM THE NB-NA AND NB-NC INTERATOMIC BOND VECTORS WITH P. C PXBA = P(1,1)*XBA + P(1,2)*YBA + P(1,3)*ZBA PYBA = P(2,1)*XBA + P(2,2)*YBA + P(2,3)*ZBA PZBA = P(3,1)*XBA + P(3,2)*YBA + P(3,3)*ZBA C PXBC = P(1,1)*XBC + P(1,2)*YBC + P(1,3)*ZBC PYBC = P(2,1)*XBC + P(2,2)*YBC + P(2,3)*ZBC PZBC = P(3,1)*XBC + P(3,2)*YBC + P(3,3)*ZBC ELSE C C NO HOUSEHOLDER TRANSFORMATION IS NECESSARY. C ITRAN = 0 PXBA = XBA PYBA = YBA PZBA = ZBA PXBC = XBC PYBC = YBC PZBC = ZBC ENDIF C C WHILE LOOKING DOWN THE X AXIS, BACK TOWARD THE ORIGIN, DETERMINE C THE DIHEDRAL ANGLE THETA THAT THE TRANSFORMED NB-NC BOND VECTOR C MAKES WITH THE POSITIVE Z AXIS C COSINE = PZBC/DSQRT(PYBC**2 + PZBC**2) IF(ABS(COSINE).GT.1.0D0) COSINE = SIGN(1.0D0,COSINE) THETA = (ACOS(COSINE)*DEGREE)*SIGN(1.0D0,PYBC) C C NOW DETERMINE THE DIHEDRAL ANGLE PHI THAT THE NA-I BOND VECTOR C WOULD MAKE WITH THE POSITIVE Z AXIS IN THIS ROTATED SYSTEM. C ACCOUNT FOR THE FACT THAT THE HOUSEHOLDER TRANSFORMATION CAUSES C A MIRROR-IMAGE EFFECT. C DIHED = ZMAT1(3,I) DIHDEG = DIHED*DEGREE IF(ITRAN.EQ.0)THEN PHI = THETA - DIHDEG ELSE PHI = THETA + DIHDEG ENDIF IF(PHI.GT.180.0D0) PHI = PHI - 360.0D0 IF(PHI.LT.-180.0D0) PHI = PHI + 360.0D0 C C DETERMINE THE RELATIVE Y AND Z COORDINATES OF THE NA-I BOND C VECTOR. IF NECESSARY, ASSIGN Y-Z QUADRANT: C C C (PHI=0) C C Z C | C 4 | 1 C | C (PHI=-90) -----|-----Y (PHI=90) C | C 3 | 2 C | C C (ABS(PHI)=180) C C DIFF = ABS(ABS(PHI)-90.0D0) IF(DIFF.LT.1.0D-3)THEN PZAI = 0.0D0 PYAI = SIGN(1.0D0,PHI) ELSE IF(PHI.GE.0.0D0.AND.PHI.LT.90.0D0)THEN IQUAD = 1 ELSEIF(PHI.GT.90.0D0)THEN IQUAD = 2 ELSEIF(PHI.LT.0.0D0.AND.PHI.GT.-90.0D0)THEN IQUAD = 4 ELSE IQUAD = 3 ENDIF C C USE THE QUADRANT NUMBER AND THE TANGENT FUNCTION TO GET C RELATIVE Y AND Z COORDINATES. AVOID PHI=90,-90. C IF(IQUAD.EQ.4.OR.IQUAD.EQ.3)THEN IF(PHI.LT.-85.0D0.AND.PHI.GT.-95.0D0)THEN PYAI = -1.0D0 PHI90 = PHI + 90.0D0 PZAI = TAN(PHI90*RADIAN) ELSEIF(IQUAD.EQ.4)THEN PZAI = 1.0D0 PYAI = -TAN(ABS(PHI*RADIAN)) ELSEIF(IQUAD.EQ.3)THEN PZAI = -1.0D0 PHI90 = 180.0D0 - ABS(PHI) PYAI = -TAN(PHI90*RADIAN) ENDIF ELSEIF(IQUAD.EQ.1.OR.IQUAD.EQ.2)THEN IF(PHI.GT.85.0D0.AND.PHI.LT.95.0D0)THEN PYAI = 1.0D0 PHI90 = 90.0D0 - PHI PZAI = TAN(PHI90*RADIAN) ELSEIF(IQUAD.EQ.1)THEN PZAI = 1.0D0 PYAI = TAN(PHI*RADIAN) ELSEIF(IQUAD.EQ.2)THEN PZAI = -1.0D0 PHI180 = 180.0D0 - PHI PYAI = TAN(PHI180*RADIAN) ENDIF ENDIF ENDIF C C ASSUMING A UNIT Y-Z COMPONENT, COMPUTE THE CORRESPONDING C X COMPONENT OF NA-I USING THE CURRENT BOND ANGLE (ALL IN C THE ROTATED SYSTEM). C YZNORM = DSQRT(PYAI**2 + PZAI**2) PYAI = PYAI/YZNORM PZAI = PZAI/YZNORM ANG90 = ANGDEG - 90.0D0 PXAI = TAN(ANG90*RADIAN) C C NOW NORMALIZE THE ROTATED NA-I VECTOR WITH RESPECT TO AN ORIGIN C AT ATOM NA. C RAI = ZMAT1(1,I) XSCALE = RAI/DSQRT(PXAI**2 + 1.0D0) PXAI = PXAI*XSCALE + PXBA PYAI = PYAI*XSCALE PZAI = PZAI*XSCALE C C TRANSFORM BACK TO THE ORIGINAL COORDINATE SYSTEM -- THE HOUSEHOLDER C MATRIX IS ITS OWN INVERSE. C IF(ITRAN.EQ.1)THEN XAI = P(1,1)*PXAI + P(1,2)*PYAI + P(1,3)*PZAI YAI = P(2,1)*PXAI + P(2,2)*PYAI + P(2,3)*PZAI ZAI = P(3,1)*PXAI + P(3,2)*PYAI + P(3,3)*PZAI ELSE XAI = PXAI YAI = PYAI ZAI = PZAI ENDIF C C CORRECT FOR THE FACT THAT THE ORIGIN FOR THE HOUSEHOLDER C TRANSFORMATION WAS ATOM ATOM B. C XYZ1(1,I) = XAI + XYZ1(1,NB) XYZ1(2,I) = YAI + XYZ1(2,NB) XYZ1(3,I) = ZAI + XYZ1(3,NB) 500 CONTINUE 1000 RETURN END