SUBROUTINE ORIENT( GEOM1, IE1, GEOM2, IE2, DISPL, ANGLES) * * COMPUTE ANGLE AND DISPLACEMENT RELATION BETWEEN GEOM1 AND GEOM2 * THE RESULT WILL BE DISPL AND ANGLES WHICH WILL CONVERT * GEOM1 INTO GEOM2 * * THE ASSUMPTION IS THAT GEOM1 IS NOT ROTATED * GEOM1( *, 1): 0.000 0.000 0.000 * GEOM1( *, 2): X1 0.000 0.000 * GEOM1( *, 3): X2 Y2 0.000 * GEOM1( *, 4): X3 Y3 Z3 * * The method here is to step-by-step find the displacements and * rotation angles to return GEOM2 back to the starting orientation * represented by GEOM1. * Throughout the part the rotation routines will be called in such * a manner as to pretend the entire GEOM2 array (here used as GEOML * a logcal copy) is a sub-group of a larger geometry. In this way, * the TOTROT information will not be accumulated. * At the end, the process will be reversed and the rotation * information will be accumulatd via TOTROT. * IMPLICIT REAL( A-H, O-Z) INCLUDE 'SIZES' DIMENSION IE1( NUMATM), IE2( NUMATM) DIMENSION GEOM1( 3,NUMATM), GEOM2( 3,NUMATM), GEOML( 3, NUMATM) DIMENSION DISPL( 3), ANGLES( 3), IGROUP( NUMATM) LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI * IF ( DEBUG) THEN WRITE (*,*) ' ORIENT: geom1 geom2' DO 9 I=1, 5 WRITE (*, '( I3, I4, 3F10.5, 4X, I4, 3F10.5)') I, IE1(I), . (GEOM1( K, I),K=1,3), IE2( I), (GEOM2( K, I),K=1,3) 9 CONTINUE ENDIF NTOTAL = 0 DO 1 I=1, NUMATM IF ( IE1(I).NE.0 .OR. IE2(I).NE.0) THEN IGROUP( I) = I NTOTAL = NTOTAL + 1 ELSE GOTO 2 ENDIF 1 CONTINUE * 2 CONTINUE I1 = 1 I2 = 1 10 CONTINUE IF( IE1( I1) .GE. 99 ) THEN I1 = I1 + 1 GOTO 10 ENDIF IF( IE2( I2) .GE. 99 ) THEN I2 = I2 + 1 GOTO 10 ENDIF IF ( IE1( I1) .NE. IE2( I2) ) THEN CALL DEBUGR( 'ORIENT: FIRST ATOMS DO NOT MATCH.') RETURN ENDIF DO 100 J= 1, 3 DISPL( J) = GEOM1( J, I1) - GEOM2( J, I2) 100 CONTINUE DO 105 K= 1, NTOTAL DO 105 J= 1, 3 GEOML( J, K) = GEOM2( J, K) + DISPL( J) 105 CONTINUE J1 = I1+1 J2 = I2+1 110 CONTINUE IF( IE1( J1) .GE. 99 ) THEN J1 = J1 + 1 GOTO 110 ENDIF IF( IE2( J2) .GE. 99 ) THEN J2 = J2 + 1 GOTO 110 ENDIF IF ( IE1( J1) .NE. IE2( J2) ) THEN CALL DEBUGR( 'ORIENT: SECOND ATOMS DO NOT MATCH.') RETURN ENDIF PI = 2.0D0 * ASIN(1.0D0) CONV = PI / 180.0D0 * FIRST ROTATE ABOUT Z AXIS TO ACHIEVE VERTICAL COINCIDENCE D12 = ( GEOM1( 1, I1) - GEOM1( 1, J1))**2 + . ( GEOM1( 2, I1) - GEOM1( 2, J1))**2 D13 = ( GEOM1( 1, I1) - (GEOML( 1, J2)))**2 + . ( GEOM1( 2, I1) - (GEOML( 2, J2)))**2 D23 = ( GEOM1( 1, J1) - (GEOML( 1, J2)))**2 + . ( GEOM1( 2, J1) - (GEOML( 2, J2)))**2 XY = SQRT( D12*D13) IF( XY.LT.1.0D-2) THEN CALL DEBUGR( 'ORIENT: ATOMS HAVE ZERO SPACING.' ) RETURN ENDIF TEMP = 0.5D0 * (D12+D13-D23) / XY IF(TEMP.GT.1.0D0) THEN TEMP=1.0D0 ELSEIF(TEMP.LT.-1.0D0)THEN TEMP=-1.0D0 ENDIF C?? ANGLES( 1) = -ACOS( TEMP) / CONV ANGLES( 1) = ACOS( TEMP) / CONV * USE CARTESIAN ROTATION ROUTINE WITH NGROUP SET CALL RCART( NTOTAL, GEOML, GEOML, I2, 0.0D0, 0.0D0, ANGLES(1), . NTOTAL, IGROUP) * NOW WORK-OUT THE ROTATION ABOUT Y-AXIS TO MAKE THEM MATCH D23 = (GEOM1( 3, I2) - GEOML( 3, J2))**2 TEMP = 0.5D0 * (D12+D13-D23) / XY IF(TEMP.GT.1.0D0) THEN TEMP=1.0D0 ELSEIF(TEMP.LT.-1.0D0)THEN TEMP=-1.0D0 ENDIF ANGLES( 2) = -ACOS( TEMP) / CONV CALL RCART( NTOTAL, GEOML, GEOML, I2, 0.0D0, ANGLES( 2), 0.0D0, . NTOTAL, IGROUP) 200 CONTINUE K1 = J1 + 1 K2 = J2 + 1 210 CONTINUE IF( IE1( K1) .GE. 99 ) THEN K1 = K1 + 1 GOTO 210 ENDIF IF( IE2( K2) .GE. 99 ) THEN K2 = K2 + 1 GOTO 210 ENDIF IF ( IE1( K1) .NE. IE2( K2) ) THEN CALL DEBUGR( 'ORIENT: SECOND ATOMS DO NOT MATCH.') RETURN ENDIF * THE FOLLOWING CODE IS MODIFIED FROM ROUTINE DIHED * ATOM I IS GEOML( X, K2) * J GEOM1( X, J1) * K GEOM1( X, I1) * L GEOM1( X, K1) XI1 = GEOML( 1, K2) - GEOM1( 1, I1) XJ1 = GEOM1( 1, J1) - GEOM1( 1, I1) XL1 = GEOM1( 1, K1) - GEOM1( 1, I1) YI1 = GEOML( 2, K2) - GEOM1( 2, I1) YJ1 = GEOM1( 2, J1) - GEOM1( 2, I1) YL1 = GEOM1( 2, K1) - GEOM1( 2, I1) ZI1 = GEOML( 3, K2) - GEOM1( 3, I1) ZJ1 = GEOM1( 3, J1) - GEOM1( 3, I1) ZL1 = GEOM1( 3, K1) - GEOM1( 3, I1) DIST = SQRT(XJ1**2+YJ1**2+ZJ1**2) COSA = ZJ1 / DIST IF(COSA.GT.1.0D0) COSA = 1.0D0 IF(COSA.LT.-1.0D0) COSA = -1.0D0 DDD=1.0D0-COSA**2 IF(DDD.LE.0.0) GOTO 315 YXDIST=DIST* SQRT(DDD) IF(YXDIST.GT.1.0D-9) GOTO 310 315 CONTINUE XI2=XI1 XL2=XL1 YI2=YI1 YL2=YL1 COSTH=COSA SINTH=0.D0 GOTO 311 310 COSPH=YJ1/YXDIST SINPH=XJ1/YXDIST XI2=XI1*COSPH-YI1*SINPH XJ2=XJ1*COSPH-YJ1*SINPH XL2=XL1*COSPH-YL1*SINPH YI2=XI1*SINPH+YI1*COSPH YJ2=XJ1*SINPH+YJ1*COSPH YL2=XL1*SINPH+YL1*COSPH * ROTATE KJ AROUND THE X AXIS SO KJ LIES ALONG THE Z AXIS COSTH=COSA SINTH=YJ2/DIST 311 CONTINUE YI3=YI2*COSTH-ZI1*SINTH YL3=YL2*COSTH-ZL1*SINTH CALL DANG(XL2,YL3,XI2,YI3,ANGLES(3)) c? ANGLES(3)=-ANGLES(3) c? IF (ANGLES(3) .LT. -180.0D0) ANGLES(3) = ANGLES(3) + 360.0D0 c? IF (ANGLES(3) .GT. 180.0D0) ANGLES(3) = ANGLES(3) - 360.0D0 CALL RPAIR( NTOTAL, GEOML, GEOML, I2, J2, ANGLES(3), . NTOTAL, IGROUP) * Now that all components are known, do the same thing in reverse and * allow the rotation routines to accumulate the info in TOTROT. DO 400 J= 1, 3 ANGLES( J) = -ANGLES( J) 400 CONTINUE DO 410 K= 1, NTOTAL DO 410 J= 1, 3 GEOML( J, K) = GEOM1( J, K) 410 CONTINUE CALL RPAIR( NTOTAL, GEOML, GEOML, I1, J1, ANGLES(3), . 0, IGROUP) CALL RCART( NTOTAL, GEOML, GEOML, I1, 0.0D0, ANGLES(2), 0.0D0, . 0, IGROUP) CALL RCART( NTOTAL, GEOML, GEOML, I1, 0.0D0, 0.0D0, ANGLES(1), . 0, IGROUP) DO 440 K= 1, NTOTAL DO 440 J= 1, 3 GEOML( J, K) = GEOML( J, K) - DISPL( J) 440 CONTINUE c? WRITE (*,*) ' ORIENT: geom1 geom2' c? M1 = 0 c? M2 = 0 c? DO 499 I=1, 5 c? M1 = M1 + 1 c? M2 = M2 + 1 c? 491 IF ( IE1( M1).EQ.99) THEN c? M1 = M1 + 1 c? GOTO 491 c? ENDIF c? 492 IF ( IE2( M2).EQ.99) THEN c? M2 = M2 + 1 c? GOTO 492 c? ENDIF c? WRITE (*, '( I3, I4, 3F10.5, 4X, I4, 3F10.5)') c? . M1,IE1(M1),(GEOML(K,M1),K=1,3),IE2(M2),(GEOM2(K,M2),K=1,3) c? 499 CONTINUE 600 CONTINUE DO 610 I= 1, NUMATM DO 610 J= 1, 3 GEOM1( J, I) = GEOML( J, I) 610 CONTINUE RETURN END