root / examples / Test / PES_HCN_ana.f90 @ 12
Historique | Voir | Annoter | Télécharger (3,99 ko)
| 1 | 1 | pfleura2 | PROGRAM PES2D |
|---|---|---|---|
| 2 | 1 | pfleura2 | |
| 3 | 1 | pfleura2 | IMPLICIT NONE |
| 4 | 1 | pfleura2 | integer, parameter :: KINT = kind(1) |
| 5 | 1 | pfleura2 | integer, parameter :: KREAL = kind(1.0d0) |
| 6 | 1 | pfleura2 | |
| 7 | 1 | pfleura2 | REAL(KREAL) :: E |
| 8 | 1 | pfleura2 | |
| 9 | 1 | pfleura2 | ! Bohr --> Angstr |
| 10 | 1 | pfleura2 | real(KREAL), parameter :: BOHR = 0.52917726D+00 |
| 11 | 1 | pfleura2 | ! |
| 12 | 1 | pfleura2 | ! Parameters to define the surface |
| 13 | 1 | pfleura2 | INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/) |
| 14 | 1 | pfleura2 | INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/) |
| 15 | 1 | pfleura2 | REAL(KREAL), PARAMETER :: hh=0.01d0 |
| 16 | 1 | pfleura2 | |
| 17 | 1 | pfleura2 | ! Variables |
| 18 | 1 | pfleura2 | INTEGER(KINT) :: i,j,iat,jat |
| 19 | 1 | pfleura2 | REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:) |
| 20 | 1 | pfleura2 | REAL(KREAL) :: xp,yp |
| 21 | 1 | pfleura2 | |
| 22 | 1 | pfleura2 | CHARACTER(132) :: Line |
| 23 | 1 | pfleura2 | LOGICAL :: Debug |
| 24 | 1 | pfleura2 | |
| 25 | 1 | pfleura2 | real(KREAL), parameter :: zero = 0.0_KREAL, one = 1.0_KREAL |
| 26 | 1 | pfleura2 | REAL(KREAL), external :: Eint |
| 27 | 1 | pfleura2 | |
| 28 | 1 | pfleura2 | REAL(KREAL) :: dCH,dCN,aHCN |
| 29 | 1 | pfleura2 | REAL(KREAL), PARAMETER :: dCHmin=1., dCHmax=2.5 |
| 30 | 1 | pfleura2 | REAL(KREAL), PARAMETER :: aHCNmin=0., aHCNmax=180. |
| 31 | 1 | pfleura2 | REAL(KREAL) :: StepdCH, StepaHCN |
| 32 | 1 | pfleura2 | INTEGER(KINT) :: NbStpdCH,NbStpaHCN |
| 33 | 1 | pfleura2 | |
| 34 | 1 | pfleura2 | |
| 35 | 1 | pfleura2 | |
| 36 | 1 | pfleura2 | ! For now dCN=1.16 is fixed. Next version should optimize it at each step ! |
| 37 | 1 | pfleura2 | dCN=1.16 |
| 38 | 1 | pfleura2 | NbStpdCH=100 |
| 39 | 1 | pfleura2 | NbStpaHCN=100 |
| 40 | 1 | pfleura2 | StepdCH=(dCHmax-dCHmin)/(NbStpdCH-1) |
| 41 | 1 | pfleura2 | StepaHCN=(aHCNmax-aHCNmin)/(NbStpaHCN-1) |
| 42 | 1 | pfleura2 | dCH=dCHmin |
| 43 | 1 | pfleura2 | DO I=1,NbStpdCH |
| 44 | 1 | pfleura2 | aHCN=aHCNmin |
| 45 | 1 | pfleura2 | DO J=1,NbStpaHCN |
| 46 | 1 | pfleura2 | e=Eint(dCH,dCN,aHCN) |
| 47 | 1 | pfleura2 | WRITE(*,'(1X,F10.4,1X,F10.4,1X,F8.2,1X,F15.6)') dCH,dCN,aHCN,e |
| 48 | 1 | pfleura2 | aHCN=aHCN+StepaHCN |
| 49 | 1 | pfleura2 | END DO |
| 50 | 1 | pfleura2 | WRITE(*,*) |
| 51 | 1 | pfleura2 | dCH=dCH+StepdCH |
| 52 | 1 | pfleura2 | END DO |
| 53 | 1 | pfleura2 | |
| 54 | 1 | pfleura2 | !==================================================================== |
| 55 | 1 | pfleura2 | end |
| 56 | 1 | pfleura2 | |
| 57 | 1 | pfleura2 | |
| 58 | 1 | pfleura2 | function Eint(dCH,dCN,aHCN) |
| 59 | 1 | pfleura2 | |
| 60 | 1 | pfleura2 | IMPLICIT NONE |
| 61 | 1 | pfleura2 | integer, parameter :: KINT = kind(1) |
| 62 | 1 | pfleura2 | integer, parameter :: KREAL = kind(1.0d0) |
| 63 | 1 | pfleura2 | |
| 64 | 1 | pfleura2 | |
| 65 | 1 | pfleura2 | REAL(KREAL) ,INTENT(IN) :: dCH,dCN,aHCN |
| 66 | 1 | pfleura2 | REAL(KREAL) :: Eint |
| 67 | 1 | pfleura2 | |
| 68 | 1 | pfleura2 | INTEGER(KINT) :: i |
| 69 | 1 | pfleura2 | REAL(KREAL) :: dNH |
| 70 | 1 | pfleura2 | |
| 71 | 1 | pfleura2 | real(KREAL) :: autoA,autoeV,Pi |
| 72 | 1 | pfleura2 | parameter (autoA=0.52917715d0) |
| 73 | 1 | pfleura2 | parameter (autoeV=27.21183d0) |
| 74 | 1 | pfleura2 | |
| 75 | 1 | pfleura2 | |
| 76 | 1 | pfleura2 | real(KREAL) :: Q(3) |
| 77 | 1 | pfleura2 | real(KREAL) :: TT,gR,pR |
| 78 | 1 | pfleura2 | |
| 79 | 1 | pfleura2 | real(KREAL) :: mB,mC,mH |
| 80 | 1 | pfleura2 | real(KREAL) :: T,RN,RC,r1,r2,r3 |
| 81 | 1 | pfleura2 | |
| 82 | 1 | pfleura2 | real(KREAL) :: Z1,Z12,V1 |
| 83 | 1 | pfleura2 | real(KREAL) :: Z2,Z22,V2 |
| 84 | 1 | pfleura2 | real(KREAL) :: Z3,Z32,V3 |
| 85 | 1 | pfleura2 | |
| 86 | 1 | pfleura2 | real(KREAL) :: S1,S2,S3,S12,S22,S32,POLY |
| 87 | 1 | pfleura2 | real(KREAL) :: E1,E2,E3,HY1,HY2,HY3,SWITCH,V123 |
| 88 | 1 | pfleura2 | |
| 89 | 1 | pfleura2 | pi=dacos(-1.0d0) |
| 90 | 1 | pfleura2 | |
| 91 | 1 | pfleura2 | dNH=sqrt(dCH**2+dCN**2-2*dCH*dCN*cos(aHCN*pi/180.)) |
| 92 | 1 | pfleura2 | |
| 93 | 1 | pfleura2 | R1=dCH |
| 94 | 1 | pfleura2 | R2=dCN |
| 95 | 1 | pfleura2 | R3=dNH |
| 96 | 1 | pfleura2 | |
| 97 | 1 | pfleura2 | |
| 98 | 1 | pfleura2 | |
| 99 | 1 | pfleura2 | ! R1=R1*autoA |
| 100 | 1 | pfleura2 | ! R2=R2*autoA |
| 101 | 1 | pfleura2 | ! R3=R3*autoA |
| 102 | 1 | pfleura2 | |
| 103 | 1 | pfleura2 | !.....CH |
| 104 | 1 | pfleura2 | Z1=R1-1.0823 |
| 105 | 1 | pfleura2 | Z12=Z1*Z1 |
| 106 | 1 | pfleura2 | V1=-2.8521*(1.0+5.5297*Z1+8.7166*Z12+5.3082*Z1*Z12)*dexp(-5.5297*Z1) |
| 107 | 1 | pfleura2 | !.....CN |
| 108 | 1 | pfleura2 | Z2=R2-1.1718 |
| 109 | 1 | pfleura2 | Z22=Z2*Z2 |
| 110 | 1 | pfleura2 | V2=-7.9282*(1.0+5.2448*Z2+7.3416*Z22+4.9785*Z2*Z22)*dexp(-5.2448*Z2) |
| 111 | 1 | pfleura2 | !.....NH |
| 112 | 1 | pfleura2 | Z3=R3-1.0370 |
| 113 | 1 | pfleura2 | Z32=Z3*Z3 |
| 114 | 1 | pfleura2 | V3=-3.9938*(1.0+3.0704*Z3)*dexp(-3.0704*Z3) |
| 115 | 1 | pfleura2 | |
| 116 | 1 | pfleura2 | ! write(6,*) v1,v2,v3 |
| 117 | 1 | pfleura2 | |
| 118 | 1 | pfleura2 | !.....THREE BODY TERMS |
| 119 | 1 | pfleura2 | Z1=R1-1.9607 |
| 120 | 1 | pfleura2 | Z2=R2-2.2794 |
| 121 | 1 | pfleura2 | Z3=R3-1.8687 |
| 122 | 1 | pfleura2 | S1=0.4436*Z1+0.6091*Z2+0.6575*Z3 |
| 123 | 1 | pfleura2 | S2=-.8941*Z1+0.2498*Z2+0.3718*Z3 |
| 124 | 1 | pfleura2 | S3=0.0622*Z1-0.7527*Z2+0.6554*Z3 |
| 125 | 1 | pfleura2 | |
| 126 | 1 | pfleura2 | S12=S1*S1 |
| 127 | 1 | pfleura2 | S22=S2*S2 |
| 128 | 1 | pfleura2 | S32=S3*S3 |
| 129 | 1 | pfleura2 | POLY=-3.0578*(1.0+1.9076*S1-0.5008*S2-0.0149*S3+0.6695*S12 & |
| 130 | 1 | pfleura2 | -1.3535*S22-1.0501*S32+0.2698*S1*S2-1.1120*S1*S3 & |
| 131 | 1 | pfleura2 | +1.9310*S2*S3-0.0877*S1*S12+0.0044*S2*S22+0.0700 & |
| 132 | 1 | pfleura2 | *S3*S32+0.0898*S12*S2-1.0186*S1*S22-0.0911*S12 & |
| 133 | 1 | pfleura2 | *S3+0.0017*S1*S32+0.4567*S22*S3-0.8840*S2*S32 & |
| 134 | 1 | pfleura2 | +0.3333*S1*S2*S3-0.0367*S12*S12+0.4821*S22*S22 & |
| 135 | 1 | pfleura2 | +0.2564*S32*S32-0.0017*S12*S1*S2-0.2278*S12*S22 & |
| 136 | 1 | pfleura2 | -0.1287*S1*S2*S22+0.1759*S1*S12*S3-0.0399*S12*S32 & |
| 137 | 1 | pfleura2 | -0.1447*S1*S3*S32-0.3147*S2*S22*S3+0.1233*S22*S32 & |
| 138 | 1 | pfleura2 | +0.3161*S2*S3*S32+0.0919*S12*S2*S3-0.0954*S1*S22*S3 & |
| 139 | 1 | pfleura2 | +0.1778*S1*S2*S32-0.1892*S22*S22*S2) |
| 140 | 1 | pfleura2 | |
| 141 | 1 | pfleura2 | |
| 142 | 1 | pfleura2 | E1=dexp(3.9742*Z1/2.) |
| 143 | 1 | pfleura2 | E2=dexp(4.3688*Z2/2.) |
| 144 | 1 | pfleura2 | E3=dexp(1.5176*Z3/2.) |
| 145 | 1 | pfleura2 | HY1=(E1-1.0/E1)/(E1+1.0/E1) |
| 146 | 1 | pfleura2 | HY2=(E2-1.0/E2)/(E2+1.0/E2) |
| 147 | 1 | pfleura2 | HY3=(E3-1.0/E3)/(E3+1.0/E3) |
| 148 | 1 | pfleura2 | SWITCH=(1.0-HY1)*(1.0-HY2)*(1.0-HY3) |
| 149 | 1 | pfleura2 | V123=SWITCH*POLY |
| 150 | 1 | pfleura2 | |
| 151 | 1 | pfleura2 | ! write(6,*) v123,V1+V2+V3+V123 |
| 152 | 1 | pfleura2 | |
| 153 | 1 | pfleura2 | Eint=(V1+V2+V3+V123)/autoeV |
| 154 | 1 | pfleura2 | |
| 155 | 1 | pfleura2 | |
| 156 | 1 | pfleura2 | |
| 157 | 1 | pfleura2 | return |
| 158 | 1 | pfleura2 | end function Eint |
| 159 | 1 | pfleura2 |