C ******************************************************************* C *************** CLASSICAL TRAJECTORY PROGRAM - 2D ************* C Copyright (C) 2002 J.M.C. Marques, A. Riganelli, A.J.C. Varandas C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the License, or C (at your option) any later version. C C This program is distributed in the hope that it will be useful, but C WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU C General Public License for more details. C C You should have received a copy of the GNU General Public License C along with this program; if not, write to the Free Software Foundation, C Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. C C ******************************************************************* C THIS PROGRAM SHOULD BE COMPILED WITH THE FOLLOWING ROUTINES: C 1) POTENTIAL OF THE SYSTEM: PES(R,V), WHERE R IS THE VECTOR C OF THE INTERNUCLEAR DISTANCES AND V IS THE POTENTIAL. THE C POTENTIAL SHOULD BE SHIFTED, SO THAT THE MINIMUM OF V FOR C THE REACTANT (A+BC) ASYMPTOTIC CHANNEL BECOMES ZERO. C 2) 1ST DERIVATIVES OF THE POTENTIAL: DPES(R,DER), WHERE R IS C THE VECTOR OF THE INTERNUCLEAR DISTANCES AND DER IS THE C VECTOR OF PARTIAL DERIVATIVES. C 3) POTENTIAL FUNCTIONS OF DIATOMICS: VDIAT(RI,V2,I), WHERE C RI IS THE INTERNUCLEAR DISTANCE OF THE Ith DIATOMIC (I=1, C 2, 3 STANDS FOR AB, BC, AC, RESPECTIVELY) AND V2 IS THE C CORRESPONDING DIATOMIC POTENTIAL. THIS MUST BE SHIFTED TO C HAVE V2=0 FOR THE CORRESPONDING EQUILIBRIUM GEOMETRY. C C ALL QUANTITIES IN THIS PROGRAM ARE IN ATOMIC UNITS. C HOWEVER, ENERGIES AND MASSES ARE READ (FROM DATA FILE) AND WRITTEN C IN KCAL/MOL AND AMU, RESPECTIVELY. C C == INPUT FILE: TRAJ2D.DAT == C DATA FILE (TRAJ2D.DAT) FOLLOWS THE ORDER: C -- MAIN PROGRAM -- C A) NUMBER OF TRAJECTORIES TO BE CALCULATED (NTRAJ). C B) INTEGRATION STEP (STEP) IN ATOMIC UNITS OF TIME; MAXIMUM C NUMBER OF TIME STEPS (NSTEP). C C) MASSES OF THE ATOMS (AMM(3)). C -- SUBROUTINE ICOND1 -- C D) EQUILIBRIUM INTERNUCLEAR DISTANCE OF THE REACTANT BC C DIATOMIC (REBC). C E) POTENTIAL ENERGY OF PRODUCTS C+AB (DELTAB) AND B+AC (DELTAC) C AS REFERRED TO REACTANTS. C F) INITIAL INTERFRAGMENTS DISTANCE (RHO); MAXIMUM IMPACT C PARAMETER (BMAX). C G) INITIAL TRANSLATIONAL ENERGY (ETR); INITIAL VIBRATIONAL C ENERGY (EVIB); INITIAL ROTATIONAL ENERGY (EROT). C C == OUTPUT FILE: TRAJ2D.RES == C THIS IS A GENERAL SELF-EXPLAINING OUTPUT FILE. IT CONTAINS MOST OF C THE INFORMATION ARISING IN THE TRAJECTORY RUN. AT THE END, IT HAS C AN HISTOGRAM SECTION WITH THE DISTRIBUTIONS OF THE TRAJECTORIES C ACCORDINGLY TO SPECIFIC VARIABLES (IMPACT PARAMETER, SCATTERING C ANGLE, ETC.) FOR EACH CHANNEL; ALSO THE CROSS SECTIONS AND THE C RESPECTIVE STANDARD DEVIATIONS ARE WRITTEN. C C THE FOLLOWING OUTPUT FILES APPEAR ONLY IN A SINGLE TRAJECTORY C CALCULATION (NTRAJ=1): C == OUTPUT FILE: TRAJ2D_COORD.RES == C JACOBI ATOM-DIATOM COORDINATES (COORDINATES OF THE PROGRAM) VS C TIME (ALL IN ATOMIC UNITS). C C == OUTPUT FILE: TRAJ2D_MOM.RES == C CONJUGATE MOMENTA OF THE JACOBI ATOM-DIATOM COORDINATES VS C TIME (ALL IN ATOMIC UNITS). C C == OUTPUT FILE: TRAJ2D_DIST.RES == C INTERNUCLEAR DISTANCES VS TIME (ALL IN ATOMIC UNITS). C C == OUTPUT FILE: TRAJ2D_ENER.RES == C KINETIC ENERGY, POTENTIAL ENERGY, TOTAL ENERGY VS TIME AND NUMBER C OF STEPS (ALL IN ATOMIC UNITS). C ******************************************************************* PROGRAM TRAJ2D IMPLICIT REAL*8 (A-H,O-Z) C SET THE SYSTEM PARAMETERS: NUMBER OF COLLIDING ATOMS (NAT), C NUMBER OF COORDINATES (NCOORD), NUMBER OF HAMILTON EQUATIONS C OF MOTION (NHAM). PARAMETER (NAT=3,NCOORD=4,NHAM=8) !NCOORD=2*NAT-2,NHAM=2*NCOORD DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) DIMENSION XI(NCOORD),PXI(NCOORD) DIMENSION IHB(0:3,20),IHTH(0:3,20),IHANG(0:3,20),ICH(0:3), & TOTCS(0:3),STD(0:3),IHEE(0:3,20),IHEV(0:3,20), & IHER(0:3,20),ETH(3) COMMON/TEMP/NSTEP COMMON/TEMP1/STEP COMMON/B13/B,THETA COMMON/B14/NTRAJ COMMON/EKEEP/ETOTI,DEMAX DATA IRAND/1361589/ OPEN(10,FILE='TRAJ2D.DAT',STATUS='OLD') OPEN(11,FILE='TRAJ2D.RES',STATUS='UNKNOWN') OPEN(12,FILE='TRAJ2D_COORD.RES',STATUS='UNKNOWN') OPEN(13,FILE='TRAJ2D_MOM.RES',STATUS='UNKNOWN') OPEN(14,FILE='TRAJ2D_DIST.RES',STATUS='UNKNOWN') OPEN(15,FILE='TRAJ2D_ENER.RES',STATUS='UNKNOWN') READ(10,*)NTRAJ READ(10,*)STEP,NSTEP READ(10,*)(AMM(I),I=1,NAT) C TRANSFORMATION OF UNITS OF MASS FROM AMU INTO ATOMIC UNITS DO I=1,NAT AMM(I)=AMM(I)*1822.8885D0 ENDDO C CALCULATING INITIAL CONDITIONS (COORDINATES AND MOMENTA): C ICOND1->FOR THE WHOLE SET OF TRAJECTORIES C ICOND2->FOR EACH PARTICULAR TRAJECTORY CALL ICOND1(X,PX,AMM,RHO,IRAND,NAT,NCOORD) WRITE(11,100) 100 FORMAT(94X,'CONSERVATION',/, & ' TRAJECT.',2X,'CHANNEL',2X,'ATTACK ANGLE',3X,'EVIB',6X, & 'EROT',6X,'EINT',4X,'ETRANS',6X, & 'EPOT',2X,'SCAT. ANGLE',2X,'ETOTF-ETOTI') C SETTING ARRAYS FOR HISTOGRAMS AND FINAL DATA ANALYSIS TO ZERO IC=0 DO I=0,3 ICH(I)=0 TOTCS(I)=0.0D0 STD(I)=0.0D0 DO J=1,20 IHB(I,J)=0 IHTH(I,J)=0 IHANG(I,J)=0 IHEE(I,J)=0 IHEV(I,J)=0 IHER(I,J)=0 ENDDO ENDDO C BEGIN THE LOOP FOR TRAJECTORIES DO K=1,NTRAJ write(*,*)'Traj= ',K,' random seed = ',irand CALL ICOND2(X,PX,AMM,RHO,IRAND,K,NTRAJ,NAT,NCOORD) C BOOK INITIAL COORDINATES AND MOMENTA FOR FUTURE USE DO I=1,NCOORD XI(I)=X(I) PXI(I)=PX(I) ENDDO T=0.0D0 !INITIATE INTEGRATION TIME (T) C:: CALL THE ROUTINE THAT CONTROLS THE INTEGRATION OF THE TRAJECTORY CALL TRAJ(X,PX,T,STEP,RHO,AMM,NTRAJ,NSTEP,NAT,NCOORD,IC) C ANALYSE OF THE PRODUCTS CALL PRODUCTS(K,IRAND,IC,AMM,XI,PXI,X,PX,NAT,NCOORD,NHAM, & SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER,TOTCS,STD) ENDDO !End loop on trajectories C FOR NTRAJ=1 WRITES THE MAXIMUM DEVIATION FROM INITIAL TOTAL ENERGY IF(NTRAJ.EQ.1)THEN WRITE(11,*) WRITE(11,*)'ENERGY CONSERVATION:' WRITE(11,*)'MAXIMUM DEVIATION FROM THE INITIAL TOTAL ENERGY=' & ,DEMAX,' kcal mol-1' ENDIF C CALL THE ROUTINE FOR WRITING THE OUTPUT CALL WRITEHIST(SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER,TOTCS,STD) END C ******************************************************************* C ********************** INITIAL CONDITIONS ********************* C ******************************************************************* C ICOND1 ==> CALCULATION OF THE TURNING-POINTS, ROTATIONAL ANGULAR C MOMENTUM AND HALF-PERIOD OF VIBRATION FOR THE BATCH OF TRAJECTORIES C ******************************************************************* SUBROUTINE ICOND1(X,PX,AMM,RHO,IRAND,NAT,NCOORD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) COMMON/TPOINT/TURN(2) COMMON/COND/BMAX,ETR,TAUV,OMEG(2),ETH(3) COMMON/ANGM/RL COMMON/REQUIL/REBC COMMON/EKEEP/ETOTI,DEMAX DATA XACC/0.000001D0/ ! accuracy for bissection routine READ(10,*)REBC !equilibrium distance READ(10,*)DELTAB,DELTAC !delta between products and reactants energies READ(10,*)RHO,BMAX READ(10,*)ETR,EVIB,EROT C CALCULATING TOTAL ENERGY (ETOTI) FOR TRANSLATIONAL ENERGY BOXING ETOTI=ETR+EVIB+EROT C SUBTRACTING ERGICITY FROM TOTAL ENERGY: EHT(I) IS THE MAXIMUM ENERGY C CONTENT AVAILABLE FOR EACH CHANNEL. ETH(1)=ETOTI ETH(2)=ETOTI-DELTAB ETH(3)=ETOTI-DELTAC C ETH(1)=FLOAT(INT(ETH(1))+1) ETH(2)=FLOAT(INT(ETH(2))+1) ETH(3)=FLOAT(INT(ETH(3))+1) C TRANSFORMATION OF ENERGIES INTO ATOMIC UNITS. EINT IS THE INTERNAL ENERGY ETR=ETR*0.00159360124D0 EVIB=EVIB*0.00159360124D0 EROT=EROT*0.00159360124D0 EINT=EVIB+EROT ETOT=ETOTI*0.00159360124D0 DEAB=DEAB*0.00159360124D0 DEBC=DEBC*0.00159360124D0 DEAC=DEAC*0.00159360124D0 ETH(1)=ETH(1)*0.00159360124D0 ETH(2)=ETH(2)*0.00159360124D0 ETH(3)=ETH(3)*0.00159360124D0 C CALCULATION OF THE REDUCED MASS OF THE REACTANT DIATOMIC BC REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) C CALCULATION OF THE ROTATIONAL ANGULAR MOMENTUM (RL) USING THE C EQUILIBRIUM GEOMETRY OF THE REACTANT DIATOMIC RL=SQRT(2.0D0*EROT*REDBC*REBC**2) C CALCULATION OF THE TURNING-POINTS DO 5 ITURN=1,2 IF(ITURN.EQ.1)THEN X1=0.01D0 X2=REBC ELSE IF(ITURN.EQ.2)THEN X1=REBC X2=6.0D0 ENDIF C TURNING-POINT ITURN USING BISSECTION METHOD TURN(ITURN)=BISSEC(X1,X2,REDBC,EINT,XACC) C ANGULAR VELOCITY OF THE DIATOMIC IN THE TURNING-POINT ITURN OMEG(ITURN)=RL/(REDBC*TURN(ITURN)**2) 5 CONTINUE C CALCULATION OF THE HALF-PERIOD OF VIBRATION (TAUV) RTURN=TURN(1) RMAX=TURN(2) OMEGA=OMEG(1) CALL PERIODV(OMEGA,RTURN,RMAX,TAUV,AMM,NAT,NCOORD) C WRITING ON OUTPUT SOME INITIAL PARAMETERS WRITE(11,99)BMAX,ETR/0.00159360124D0,EVIB/0.00159360124D0, & EROT/0.00159360124D0,ETOT/0.00159360124D0 99 FORMAT('************* FIXED INITIAL CONDITIONS **************',/, & ' MAXIMUM IMPACT PARAMETER/a0=',F10.5,/, & ' ETRANS/kcal mol-1= ',F10.5,/, & ' EVIB/kcal mol-1= ',F10.5,/, & ' EROT/kcal mol-1= ',F10.5,/, & ' ETOT/kcal mol-1= ',F10.5,/, & '*****************************************************') RETURN END C ******************************************************************* C ICOND2 ==> CALCULATION OF THE COORDINATES AND CONJUGATE MOMENTA FOR C EACH TRAJECTORY C ******************************************************************* SUBROUTINE ICOND2(X,PX,AMM,RHO,IRAND,ITRAJ,NTRAJ,NAT,NCOORD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) COMMON/B13/B,THETA COMMON/TPOINT/TURN(2) COMMON/COND/BMAX,ETR,TAUV,OMEG(2),ETH(3) PI=ACOS(-1.0D0) REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) REDABC=AMM(1)*(AMM(2)+AMM(3))/(AMM(1)+AMM(2)+AMM(3)) IF(MOD(ITRAJ,2).EQ.1)THEN !considering iturn=1 for the 1st trajectory RTURN=TURN(1) OMEGA=OMEG(1) ELSE RTURN=TURN(2) OMEGA=OMEG(2) ENDIF RAND=RANDOM(IRAND) C IF(ITRAJ.EQ.NTRAJ) write(50,*)irand, rand THETA=2.0D0*PI*RAND X(1)=RTURN*COS(THETA) X(2)=RTURN*SIN(THETA) RAND=RANDOM(IRAND) B=SQRT(RAND*BMAX**2) X(3)=B RAND=RANDOM(IRAND) X(4)=SQRT(RHO**2-B**2)-RAND*TAUV*SQRT(2.0D0*ETR/REDABC) PX(1)=-REDBC*OMEGA*RTURN*SIN(THETA) PX(2)=REDBC*OMEGA*RTURN*COS(THETA) PX(3)=0.0D0 PX(4)=-SQRT(2.0D0*ETR*REDABC) RETURN END C ******************************************************************* C BISSEC ==> BISSECTIONS ALGORITHM TO CALCULATE THE REQUIRED TURNING- C POINT FOR A GIVEN INTERNAL ENERGY C ******************************************************************* FUNCTION BISSEC(X1,X2,REDBC,EINT,XACC) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (JMAX=40) fmid=func(x2,REDBC,EINT) f=func(x1,REDBC,EINT) if(f*fmid.ge.0.d0) pause 'root must be bracketed in bissec' if(f.lt.0.d0)then BISSEC=x1 dx=x2-x1 else BISSEC=x2 dx=x1-x2 endif do 11 j=1,JMAX dx=dx*0.5D0 xmid=BISSEC+dx fmid=func(xmid,REDBC,EINT) if(fmid.le.0.0D0)BISSEC=xmid if(abs(dx).lt.xacc .or. fmid.eq.0.d0)then return endif 11 CONTINUE PAUSE 'CONVERGENCE NOT REACHED IN 40 ITERATIONS' END C ******************************************************************* C FUNC ==> ESTABLISHES THE BISSECTIONS' FUNCTION FOR THE CALCULATION C OF THE TURNING-POINTS ASSOCIATED TO A GIVEN INTERNAL ENERGY OF THE C BC DIATOMIC C ******************************************************************* FUNCTION FUNC(R,REDM,EINT) IMPLICIT REAL*8 (A-H,O-Z) COMMON/ANGM/RL C ADDITION OF THE CENTRIFUGAL DISTORTION TO THE DIATOMIC POTENTIAL VCENT=CENTRIF(RL,REDM,R) C FOR DIATOMIC BC I=2 I=2 CALL VDIAT(R,V2,I) FUNC=VCENT+V2-EINT RETURN END C ******************************************************************* C CENTRIF ==> CALCULATES THE CENTRIFUGAL POTENTIAL ASSOCIATED WITH C THE ROTATIONAL ANGULAR MOMENTUM (RL) C ******************************************************************* FUNCTION CENTRIF(RL,REDBC,R) IMPLICIT REAL*8 (A-H,O-Z) C CALCULATION OF THE CENTRIFUGAL POTENTIAL OF THE BC DIATOMIC CENTRIF=RL**2/(2.0D0*REDBC*R**2) RETURN END C ******************************************************************* C PERIODV ==> CALCULATES THE HALPH-PERIOD OF VIBRATION (TAUV) FOR THE C REACTANT DIATOMIC (BC) BY INTEGRATION THE CORRESPONDING HAMILTON C EQUATIONS FOLLOWED BY A THREE-POINT INTERPOLATION C ******************************************************************* SUBROUTINE PERIODV(OMEGA,RTURN,RMAX,TAUV,AMM,NAT,NCOORD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(4),PX(4),AMM(NAT) DIMENSION RANT(2),TIME(2) DIMENSION R(3),Y(8) COMMON/temp/nstep COMMON/temp1/step REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) NHAM=8 !NUMBER OF HAMILTON EQUATIONS C X(1) AND X(2) DEFINE THE B-C DISTANCE; X(3) AND X(4) DEFINE THE C DISTANCE OF A TO THE CENTER OF MASS OF BC: WITHOUT LACK OF C GENERALITY, THE DIATOMIC IS INITIALLY ALONG THE X-AXIS (THETA=0.0) C AND THE THIRD ATOM IS MOTIONLESS FAR AWAY (X(4)=100.0) FROM BC. C PX(1) TO PX(4) ARE THE CORRESPONDING CONJUGATE MOMENTA. THETA=0.0D0 X(1)=RTURN X(2)=0.0D0 X(3)=0.0D0 X(4)=100.0D0 PX(1)=-REDBC*OMEGA*RTURN*SIN(THETA) PX(2)=REDBC*OMEGA*RTURN*COS(THETA) PX(3)=0.0D0 PX(4)=0.0D0 C FIXING INITIAL VALUES OF THE INTERATOMIC DISTANCE AND TIME RANT(1)=0.0D0 RANT(2)=RTURN TIME(1)=0.0D0 TIME(2)=0.0D0 T=0.0D0 DO I=1,NSTEP CALL INTEGR(X,PX,Y,T,STEP,AMM,I,NAT,NCOORD,NHAM) CALL Q2R(Y,AMM,R,NAT,NHAM) IF((I.GT.1).AND.(R(2).LT.RANT(2)))THEN C CALCULATION OF THE RELEVANT THREE-POINT INTERPOLATION PARAMATERS C (A, B) TO OBTAIN THE VALUE OF TAUV (WHICH CORRESPONDS TO THE MAXIMUM C OF THE PARABOLA A=((RANT(1)**2-RANT(2)**2)*(T-TIME(2))+(R(2)**2-RANT(2)**2) * *(TIME(2)-TIME(1)))/((TIME(2)-TIME(1))*(T-TIME(1)) * *(T-TIME(2))) B=((RANT(1)**2-RANT(2)**2)-A*(TIME(1)**2-TIME(2)**2)) * /(TIME(1)-TIME(2)) TAUV=-B/(2.0D0*A) RETURN ELSE C BOOKING THE TWO PREVIOUS VALUES OF R(2) AND INTEGRATION TIME RANT(1)=RANT(2) RANT(2)=R(2) TIME(1)=TIME(2) TIME(2)=T ENDIF ENDDO STOP 'PROBLEMS IN THE CALCULATION OF TAUV' END C ******************************************************************* C RANDOM ==> CALCULATES RANDOM NUMBERS FOR MONTE CARLO INTEGRATION C ******************************************************************* FUNCTION RANDOM(KRANDF) IMPLICIT REAL*8 (A-H,O-Z) C C THIS ALGORITHM IS EXTRACTED FROM THE A+BC TRAJECTORY PROGRAM; C QCPE 11,273(1975). IT USES A CONGRUENCE FORMULA DESIGNED TO C WORK FOR FLOATING FRACTION LENGTHS AS SHORT AS 21 BITS. IT IS C R(I+1)=R(I)*5**2N MOD2**M. ANY ODD N IS OK AND THE LARGER M C AND N THE BETTER. R(0) MUST BE ODD. C KRAND=KRANDF D=48828125.0D0 E=2097152.0D0 Z=DBLE(KRAND) A=Z*D B=DMOD(A,E) KRAND=IFIX(SNGL(B)) RANDN=B/E KRANDF=KRAND RANDOM=RANDN RETURN END C ******************************************************************* C **************** INTEGRATION OF THE TRAJECTORY **************** C ******************************************************************* C TRAJ ==> CONTROLS THE INTEGRATION OF EACH TRAJECTORY: TESTS FOR THE C END OF THE TRAJECTORY AND WRITES INTERMEDIATE RELEVANT INFORMATION C ******************************************************************* SUBROUTINE TRAJ(X,PX,T,STEP,RHO,AMM,NTRAJ,NSTEPS,NAT,NCOORD,IC) IMPLICIT REAL*8 (A-H,O-Z) COMMON/EKEEP/ETOTI,DEMAX DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) DIMENSION R(3),Y(8) C DEMAX ACCOUNTS FOR CONSERVATION OF ENERGY DURING THE INTEGRATION C OF ONE TRAJECTORY. DEMAX=0.0D0 C NHAM= usually calculated using NAT NHAM=8 !NUMBER OF HAMILTON EQUATIONS INSTEP=500 !number of initial steps to begin the test C X(1)=Y(1) AND X(2)=Y(2) DEFINE THE B-C DISTANCE; X(3)=Y(3) AND C X(4)=Y(4) DEFINE THE DISTANCE OF A TO THE CENTER OF MASS OF BC. C PX(1) TO PX(4) ARE THE CORRESPONDING CONJUGATE MOMENTA. C C WRITE HEADINGS FOR SINGLE-TRAJECTORY OUTPUT FILES IF(NTRAJ.EQ.1)THEN WRITE(12,100) WRITE(13,101) WRITE(14,102) WRITE(15,103) ENDIF 100 FORMAT(8X,'TIME',11X,'X(1)',11X,'X(2)',11X,'X(3)',11X,'X(4)') 101 FORMAT(8X,'TIME',10X,'PX(1)',10X,'PX(2)',10X,'PX(3)',10X,'PX(4)') 102 FORMAT(8X,'TIME',11X,'R(1)',11X,'R(2)',11X,'R(3)') 103 FORMAT(' STEP ', ' TIME ',10X,'EKIN',11X,'VPOT',11X,'ETOT') DO I=1,NSTEPS CALL INTEGR(X,PX,Y,T,STEP,AMM,I,NAT,NCOORD,NHAM) IF(I.GT.INSTEP)THEN CALL Q2R(Y,AMM,R,NAT,NHAM) IF((R(1).GT.RHO).AND.(R(2).GT.RHO).AND.(R(3).GT.RHO))THEN IC=3 RETURN ELSE IF((R(2).GT.RHO).AND.(R(3).GT.RHO))THEN IC=1 RETURN ELSE IF((R(1).GT.RHO).AND.(R(2).GT.RHO))THEN IC=2 RETURN ELSE IF((R(1).GT.RHO).AND.(R(3).GT.RHO))THEN IC=0 RETURN ENDIF ENDIF C IF NTRAJ=1 (FOR ONLY ONE TRAJECTORY) WRITES IN OUTPUT FILES: C COORDINATES VS TIME; C MOMENTA VS TIME; C INTERNUCLEAR DISTANCES VS TIME; C KINETIC, POTENTIAL AND TOTAL ENERGIES VS TIME. C REDBC = REDUCED MASS OF BC MOLECULE; REDABC = REDUCED MASS OF A+BC SYSTEM. IF(NTRAJ.EQ.1)THEN REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) REDABC=AMM(1)*(AMM(2)+AMM(3))/(AMM(1)+AMM(2)+AMM(3)) CALL Q2R(Y,AMM,R,NAT,NHAM) CALL PES(R,VPOT) EKIN=(PX(1)**2+PX(2)**2)/(2.0D0*REDBC)+ & (PX(3)**2+PX(4)**2)/(2.0D0*REDABC) ETOT=EKIN+VPOT C CALCULATING MAXIMUM DEVIATION FROM INITIAL TOTAL ENERGY IF(ABS(ETOTI-ETOT).GT.DEMAX)DEMAX=ABS(ETOTI-ETOT*627.509552) C:: WRITE TO CHECK ENERGY CONSERVATION EACH 100 STEPS IF(MOD(I,100).EQ.0)THEN WRITE(12,'(1X,5E15.6)')T,Y(1),Y(2),Y(3),Y(4) WRITE(13,'(1X,5E15.6)')T,PX(1),PX(2),PX(3),PX(4) WRITE(14,'(1X,4E15.6)')T,R(1),R(2),R(3) WRITE(15,'(I7,4E15.6)')I,T,EKIN*627.509552,VPOT*627.509552 & ,ETOT*627.509552 ENDIF ENDIF ENDDO WRITE(11,*)'WARNING: TRAJECTORY HAS REACHED NO END=',I-1,'STEPS' RETURN END C ******************************************************************* C INTEGR ==> USES A FOURTH-ORDER RUNGE-KUTTA METHOD TO INTEGRATE THE C HAMILTON EQUATIONS OF MOTION C ******************************************************************* SUBROUTINE INTEGR(X,PX,Y,T,STEP,AMM,ISTEP,NAT,NCOORD,NHAM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) DIMENSION C(4),A(8,4),B(2:4,3),Y(8),COORD(8),DCOORD(8) C INTEGRATION METHOD: RUNGE-KUTTA (4TH ORDER) + HAMMING PREDICTOR-CORRECTOR DATA C/0.17476028D0,-0.55148053D0,1.20553547D0,0.17118478D0/ DATA (B(2,I),I=1,3)/0.4D0,0.0D0,0.0D0/ DATA (B(3,I),I=1,3)/0.29697760D0,0.15875966D0,0.0D0/ DATA (B(4,I),I=1,3)/0.21810038D0,-3.05096470D0,3.83286432D0/ PARAMETER (IORDER=4) REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) REDABC=AMM(1)*(AMM(2)+AMM(3))/(AMM(1)+AMM(2)+AMM(3)) DO 1 I=1,NHAM/2 Y(I)=X(I) 1 CONTINUE DO 2 I=NHAM/2+1,NHAM Y(I)=PX(I-NHAM/2) 2 CONTINUE C C 4TH RUNGE-KUTTA METHOD (IORDER=4) DO 4 K=1,IORDER DO 3 J=1,NHAM COORD(J)=Y(J) DO 3 I=K-1,1,-1 COORD(J)=COORD(J)+B(K,I)*A(J,I) 3 CONTINUE IF(K.EQ.IORDER)T=T+STEP CALL HAMILT(AMM,COORD,DCOORD,NAT,NHAM) DO 4 J=1,NHAM A(J,K)=STEP*DCOORD(J) 4 CONTINUE DO 10 J=1,NHAM DO 10 I=1,IORDER Y(J)=Y(J)+C(I)*A(J,I) 10 CONTINUE do i=1,nham/2 x(i)=y(i) enddo do i=1,nham/2 px(i)=y(i+nham/2) enddo RETURN END C ******************************************************************* C HAMILT ==> ESTABLISHES THE HAMILTON EQUATIONS FOR THE SYSTEM; SINCE C THE POTENTIAL IS SUPPOSED TO DEPEND ON THE INTERNUCLEAR DISTANCES, C THE CHAIN-RULE HAS BEEN USED TO OBTAIN THE CORRESPONDING DERIVATIVES C IN ORDER TO EACH JACOBI COORDINATE C ******************************************************************* SUBROUTINE HAMILT(AMM,COORD,DCOORD,NAT,NHAM) C HAMILTON EQUATIONS FOR A LINEAR A+BC SYSTEM. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMM(NAT) DIMENSION COORD(8),DCOORD(8) DIMENSION R(3),DER(3) C REDBC = REDUCED MASS OF BC MOLECULE; REDABC = REDUCED MASS OF A+BC SYSTEM. REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) REDABC=AMM(1)*(AMM(2)+AMM(3))/(AMM(1)+AMM(2)+AMM(3)) C TIME DERIVATIVES OF COORDINATES. DCOORD(1)=COORD(5)/REDBC DCOORD(2)=COORD(6)/REDBC DCOORD(3)=COORD(7)/REDABC DCOORD(4)=COORD(8)/REDABC C TIME DERIVATIVES OF MOMENTA. DPOT CALCULATES DERIVATIVES OF THE POTENTIAL C IN ORDER TO THE INTERNAL COORDINATES (R1, R2, R3, ...). C R(1) = AB, R(2) = BC, R(3) = AC. C CENTRIF CALCULATES THE CENTRIFUGAL POTENTIAL FOR A GIVEN ROTATIONAL ANGULAR C MOMENTUM (RL): USED IN THE CALCULATION OF THE INITIAL CONDITIONS (TAUV). CALL Q2R(COORD,AMM,R,NAT,NHAM) CALL DPES(R,DER) C AUXILIAR MASS VARIABLES (RMCB AND RMBC). RMCB=AMM(3)/(AMM(2)+AMM(3)) RMBC=AMM(2)/(AMM(2)+AMM(3)) DCOORD(5)=RMCB*(COORD(3)-RMCB*COORD(1))*DER(1)/R(1) * -COORD(1)*DER(2)/R(2)-RMBC*(COORD(3)+RMBC* * COORD(1))*DER(3)/R(3) DCOORD(6)=RMCB*(COORD(4)-RMCB*COORD(2))*DER(1)/R(1) * -COORD(2)*DER(2)/R(2)-RMBC*(COORD(4)+RMBC* * COORD(2))*DER(3)/R(3) DCOORD(7)=-(COORD(3)-RMCB*COORD(1))*DER(1)/R(1)-(COORD(3) * +RMBC*COORD(1))*DER(3)/R(3) DCOORD(8)=-(COORD(4)-RMCB*COORD(2))*DER(1)/R(1)-(COORD(4) * +RMBC*COORD(2))*DER(3)/R(3) RETURN END C ******************************************************************* C Q2R ==> CALCULATION OF THE INTERNUCLEAR DISTANCES FROM THE JACOBI C COORDINATES C ******************************************************************* SUBROUTINE Q2R(COORD,AMM,R,NAT,NHAM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION COORD(8),AMM(NAT),R(3) R(1)=SQRT((COORD(3)-(AMM(3)/(AMM(2)+AMM(3)))*COORD(1))**2 * +(COORD(4)-(AMM(3)/(AMM(2)+AMM(3)))*COORD(2))**2) R(2)=SQRT(COORD(1)**2+COORD(2)**2) R(3)=SQRT((COORD(3)+(AMM(2)/(AMM(2)+AMM(3)))*COORD(1))**2 * +(COORD(4)+(AMM(2)/(AMM(2)+AMM(3)))*COORD(2))**2) RETURN END C ******************************************************************* C ********************** PRODUCTS ANALYSIS ********************** C ******************************************************************* C PRODUCTS ==> TRANSFORMATION INTO THE JACOBI COORDINATES OF THE C OUTCOME PRODUCTS FOLLOWED BY THE CALCULATION OF THE INTERNAL ENERGY C OF THE PRODUCT DIATOMIC, ATOM-DIATOM TRANSLATIONAL ENERGY, INTERNAL C ANGULAR MOMENTUM, ROTATIONAL AND VIBRATIONAL ENERGIES OF THE C DIATOMIC, COMPONENTS OF VELOCITIES, AND SCATTERING ANGLE C ******************************************************************* SUBROUTINE PRODUCTS(K,IRAND,IC,AMM,XI,PXI,X,PX,NAT,NCOORD,NHAM, & SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER, & TOTCS,STD) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 LZ COMMON/B13/B,THETA COMMON/B14/NTRAJ COMMON/EKEEP/ETOTI,DEMAX COMMON/COND/BMAX,ETR,TAUV,OMEG(2),ETH(3) DIMENSION XI(NCOORD),PXI(NCOORD) DIMENSION X(NCOORD),PX(NCOORD),AMM(NAT) DIMENSION COORD(8) DIMENSION C(3,3),FMASS(3,3),D(3,3) DIMENSION XF(4),PXF(4),R(3) DIMENSION IHB(0:3,20),IHTH(0:3,20),IHANG(0:3,20),ICH(0:3), & TOTCS(0:3),STD(0:3),IHEE(0:3,20),IHEV(0:3,20),IHER(0:3,20) REDBC=AMM(2)*AMM(3)/(AMM(2)+AMM(3)) REDABC=AMM(1)*(AMM(2)+AMM(3))/(AMM(1)+AMM(2)+AMM(3)) C BOOK RELATIVE COORDINATES AND MOMENTA AND CALCULATE FINAL C INTERNUCLEAR DISTANCES R(1), R(2), AND R(3) DO 5 I=1,NCOORD COORD(I)=X(I) COORD(I+4)=PX(I) 5 CONTINUE CALL Q2R(COORD,AMM,R,NAT,NHAM) C CALCULATE NEW SET OF ATOM-DIATOM RELATIVE COORDINATES AND MOMENTA C USING THE APPROPRIATE TRANSFORMATION MATRIX C(3x3): C=BA-1, WHERE C B(3x3) IS THE TRANSFORMATION MATRIX BETWEEN CARTESIAN COORDINATES C AND RELATIVE COORDINATES OF THE SPECIFIC REACTION CHANNEL, WHILE C A-1(3x3) IS THE INVERSE OF THE CORRESPONDING MATRIX FOR THE C REACTANTS. IF(IC.EQ.0)THEN RR=R(2) I=2 CALL VDIAT(RR,V2,I) FMASS(1,1)=1.0D0 FMASS(1,2)=0.0D0 FMASS(1,3)=0.0D0 FMASS(2,1)=0.0D0 FMASS(2,2)=1.0D0 FMASS(2,3)=0.0D0 FMASS(3,1)=0.0D0 FMASS(3,2)=0.0D0 FMASS(3,3)=1.0D0 C(1,1) = 1.0D0 C(1,2) = 0.0D0 C(1,3) = 0.0D0 C(2,1) = 0.0D0 C(2,2) = 1.0D0 C(2,3) = 0.0D0 C(3,1) = 0.0D0 C(3,2) = 0.0D0 C(3,3) = 1.0D0 ELSE IF(IC.EQ.1)THEN RR=R(1) I=1 CALL VDIAT(RR,V2,I) REDAB=AMM(2)*AMM(1)/(AMM(2)+AMM(1)) REDCAB=AMM(3)*(AMM(2)+AMM(1))/(AMM(1)+AMM(2)+AMM(3)) FMASS(1,1)=REDAB/REDBC FMASS(1,2)=REDAB/REDABC FMASS(1,3)=0.0D0 FMASS(2,1)=REDCAB/REDBC FMASS(2,2)=REDCAB/REDABC FMASS(2,3)=0.0D0 FMASS(3,1)=0.0D0 FMASS(3,2)=0.0D0 FMASS(3,3)=1.0D0 C(1,1) = -AMM(3)/(AMM(2)+AMM(3)) C(1,2) = 1.0d0 C(1,3) = 0.0D0 C(2,1) = - AMM(2)*(AMM(1)+AMM(2)+AMM(3))/(AMM(1)+AMM(2))/ * (AMM(2)+AMM(3)) C(2,2) = -AMM(1)/(AMM(1)+AMM(2)) C(2,3) = 0.0D0 C(3,1) = 0.0D0 C(3,2) = 0.0D0 C(3,3) = 1.0D0 ELSE IF(IC.EQ.2)THEN RR=R(3) I=3 CALL VDIAT(RR,V2,I) REDAC=AMM(3)*AMM(1)/(AMM(3)+AMM(1)) REDBAC=AMM(2)*(AMM(3)+AMM(1))/(AMM(1)+AMM(2)+AMM(3)) FMASS(1,1)=REDAC/REDBC FMASS(1,2)=REDAC/REDABC FMASS(1,3)=0.0D0 FMASS(2,1)=REDBAC/REDBC FMASS(2,2)=REDBAC/REDABC FMASS(2,3)=0.0D0 FMASS(3,1)=0.0D0 FMASS(3,2)=0.0D0 FMASS(3,3)=1.0D0 C(1,1) = AMM(2)/(AMM(2)+AMM(3)) C(1,2) = 1.0D0 C(1,3) = 0.0D0 C(2,1) = AMM(3)*(AMM(1)+AMM(2)+AMM(3))/(AMM(1)+AMM(3))/ * (AMM(2)+AMM(3)) C(2,2) = -AMM(1)/(AMM(1)+AMM(3)) C(2,3) = 0.0D0 C(3,1) = 0.0D0 C(3,2) = 0.0D0 C(3,3) = 1.0D0 ENDIF DO 9 I=1,3 DO 9 J=1,3 D(I,J)=C(I,J)*FMASS(I,J) 9 CONTINUE C ****** CALCULATION OF NEW COORDINATES (XF) AND MOMENTA (PXF) ****** C USING THE FIRST LINE OF C MATRIX AND D MATRIX TO CALCULATE THE C FINAL DIATOMIC JACOBI COORDINATES AND MOMENTA, RESPECTIVELY DO 10 I=1,NCOORD/2 XF(I)=0.0D0 PXF(I)=0.0D0 JJ=0 DO 10 J=I,I+2,2 JJ=JJ+1 XF(I)=XF(I)+C(1,JJ)*X(J) PXF(I)=PXF(I)+D(1,JJ)*PX(J) 10 CONTINUE C USING THE SECOND LINE OF C MATRIX AND D MATRIX TO CALCULATE THE C FINAL ATOM-DIATOMIC JACOBI COORDINATES AND MOMENTA, RESPECTIVELY DO 11 I=NCOORD/2+1,NCOORD XF(I)=0.0D0 PXF(I)=0.0D0 JJ=0 DO 11 J=I-2,I,2 JJ=JJ+1 XF(I)=XF(I)+C(2,JJ)*X(J) PXF(I)=PXF(I)+D(2,JJ)*PX(J) 11 CONTINUE C********** CALCULATION OF FINAL-STATE PROPERTIES OF PRODUCTS ************ c CALL PES(R,POT) C ************* INTERNAL ENERGY OF THE DIATOMIC FRAGMENT ************ EINT=(PXF(1)**2+PXF(2)**2)/(2.0D0*FMASS(1,1)*REDBC)+V2 C ***************** FINAL ATOM-DIATOM TRANSLATIONAL ENERGY **************** EKIN=(PXF(3)**2+PXF(4)**2)/(2.0D0*FMASS(2,2)*REDABC) C ******************** TOTAL ENERGY OF THE SYSTEM ******************* CALL PES(R,POT) ETOT=(PXF(1)**2+PXF(2)**2)/(2.0D0*FMASS(1,1)*REDBC)+ & (PXF(3)**2+PXF(4)**2)/(2.0d0*FMASS(2,2)*REDABC)+POT C ************* INTERNAL ANGULAR MOMENTUM COMPONENT Lz ************** LZ=XF(1)*PXF(2)-XF(2)*PXF(1) C ******* ROTATIONAL AND VIBRATIONAL ENERGIES OF THE DIATOMIC ******* EROT=LZ**2/(2.0D0*FMASS(1,1)*REDBC*RR**2) EVIB=EINT-EROT C ********** COMPONENTS OF INITIAL AND FINAL VELOCITIES ************* VXI=PXI(3)/REDABC VYI=PXI(4)/REDABC VXF=PXF(3)/(FMASS(2,2)*REDABC) VYF=PXF(4)/(FMASS(2,2)*REDABC) C ************************ SCATTERING ANGLE ************************* ARG=(VXI*VXF+VYI*VYF)/(SQRT(VXI**2+VYI**2)*SQRT(VXF**2+VYF**2)) SANG=ACOS(ARG) WRITE(11,101)K,IC,THETA*180.0/ACOS(-1.0D0),EVIB*627.509552, & EROT*627.509552,EINT*627.509552,EKIN*627.509552, & POT*627.509552,SANG*180.0/ACOS(-1.0D0),ETOT*627.509552-ETOTI 101 FORMAT(I5,9X,I1,6X,6(F9.4,1X),F9.4,4X,F9.4) C ******* HIST ROUTINE PUTS FINAL QUANTITIES IN HISTOGRAM FORM ******* CALL HIST(IC,SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER,TOTCS,STD,EKIN & ,EVIB,EROT,ETOT) RETURN END SUBROUTINE HIST(IC,SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER, & TOTCS,STD,EKIN,EVIB,EROT,ETOT) C************************************************************************* C*** THIS ROUTINE CALCULATES THE HISTOGRAMS *** C************************************************************************* IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NBOX=20) COMMON/B13/B,THETA COMMON/B14/NTRAJ COMMON/COND/BMAX,ETR,TAUV,OMEG(2),ETH(3) DIMENSION IHB(0:3,20),IHTH(0:3,20),IHANG(0:3,20),ICH(0:3), & TOTCS(0:3),STD(0:3),IHEE(0:3,20),IHEV(0:3,20),IHER(0:3,20) PI=ACOS(-1.0D0) BOX=NBOX BOXP=NBOX+1 DANGLE1 = 360.D0/BOX THET= (180.D0/PI)*THETA DANGLE2 = 180.D0/BOX ANGL= (180.D0/PI)*SANG ICH(IC)=ICH(IC)+1 PROB=REAL(ICH(IC))/REAL(NTRAJ) TOTCS(IC)=PI*BMAX*BMAX*PROB IF(ICH(IC).EQ.0)THEN STD(IC)=0.0D0 ELSE STD(IC)=TOTCS(IC)*SQRT((REAL(NTRAJ)-ICH(IC))/ & (REAL(NTRAJ)*ICH(IC))) ENDIF C:: BOX FOR IMPACT PARAMETER I= IDINT(B*BOX/BMAX) +1 IHB(IC,I)= IHB(IC,I) + 1 C:: BOX FOR ATTACK ANGLE (THETA) I= IDINT(THET/DANGLE1) +1 IHTH(IC,I)= IHTH(IC,I) +1 C:: BOX FOR SCATTERING ANGLE I= IDINT(ANGL/DANGLE2) +1 IHANG(IC,I)= IHANG(IC,I) +1 C:: BOX FOR TRANSLATIONAL ENERGY I= IDINT(EKIN*BOX/ETH(IC+1)) +1 IF(I .GT. BOXP) I=BOXP IHEE(IC,I)= IHEE(IC,I) +1 C:: BOX FOR VIBRATIONAL ENERGY I= IDINT(EVIB*BOX/ETH(IC+1)) +1 IF(I .GT. BOXP) I=BOXP IHEV(IC,I)= IHEV(IC,I) +1 C:: BOX FOR ROTATIONAL ENERGY I= IDINT(EROT*BOX/ETH(IC+1)) +1 IF(I .GT. BOXP) I=BOXP IHER(IC,I)= IHER(IC,I) +1 RETURN END SUBROUTINE WRITEHIST(SANG,IHB,IHTH,IHANG,ICH,IHEE,IHEV,IHER, & TOTCS,STD) C******************************************************************************* C************ THIS ROUTINE WRITES ON OUTPUT THE FINAL HISTOGRAMS *************** C******************************************************************************* IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NBOX=20) COMMON/B13/B,THETA COMMON/B14/NTRAJ COMMON/COND/BMAX,ETR,TAUV,OMEG(2),ETH(3) DIMENSION IHB(0:3,20),IHTH(0:3,20),IHANG(0:3,20),ICH(0:3), & TOTCS(0:3),STD(0:3),IHEE(0:3,20),IHEV(0:3,20),IHER(0:3,20) C WRITE(11,200) 200 FORMAT(/" FINAL HISTOGRAMS" ) DO I=0,3 !LOOP on Channels IF(I.EQ.0) THEN WRITE(11,201) ICH(I),TOTCS(I), STD(I) 201 FORMAT(/1x,'PRODUCT A+BC',5X,'NUMBER OF TJS= ',I5, & ' CROSS SECTION/a0^2=',F12.6,' STAND. DEVI.=',F12.6//) ELSE IF(I.EQ.1) THEN WRITE(11,202) ICH(I),TOTCS(I), STD(I) 202 FORMAT(//1x,'PRODUCT AB+C',5X,'NUMBER OF TJS= ',I5, & ' CROSS SECTION/a0^2=',F12.6, ' STAND. DEVI.=',F12.6/) ELSE IF(I.EQ.2) THEN WRITE(11,203) ICH(I),TOTCS(I), STD(I) 203 FORMAT(//1x,'PRODUCT AC+B',5X,'NUMBER OF TJS= ',I5, & ' CROSS SECTION/a0^2=',F12.6,' STAND. DEVI.=',F12.6/) ELSE IF(I.EQ.3) THEN WRITE(11,204) ICH(I),TOTCS(I), STD(I) 204 FORMAT(//1x,'PRODUCT A+B+C',5X,'NUMBER OF TJS= ',I5, & ' CROSS SECTION/a0^2=',F12.6,' STAND. DEVI.=',F12.6/) END IF IF(I.NE.3) THEN WRITE(11,205) ETH(I+1)*627.509552 205 FORMAT('RANGE FOR THE ENERGY HISTOGRAMS/(kcal mol-1):', & 2x,'0 - ',F8.4) WRITE(11,206) 206 FORMAT(3X,' IMP. PAR. ','ATTACK ANGLE',5x, & 'E. TRANS.',7x,'E. VIB.',8x,'E. ROT.',5x,'SCATT. ANGLE') WRITE(11,207)BMAX 207 FORMAT(' 0 - ',F5.2, ' 0 - 360',47X,' 0 - 180') DO J=1,NBOX WRITE(11,208) J,IHB(I,J),IHTH(I,J),IHEE(I,J), & IHEV(I,J),IHER(I,J),IHANG(I,J) END DO 208 FORMAT(I3,4X,I5,4X,5(4X,I5,6X)) END IF END DO RETURN END