PROGRAM
POTENTIAL
IMPLICIT
REAL*8 (A-H,O-Z)
DIMENSION
DIPPROPS(5)
C
C THIS PROGRAM
GENERATES A POTENTIAL FOR PROTON ATTACHMENT OF HCN
C AND FOR
PROTON TRANSFER IN LINEAR HCN - H+ - NCH (SCF-LEVEL MODEL).
C
C THE ENERGY EVALUATION
REQUIRES TWO DISTANCES WHICH ARE THE DISTANCES
C FROM THE
PROTON TO THE NITROGEN ATOMS IN AGNSTROMS.
C TO GENERATE
AN H+ - NCH POTENTIAL, SET THE SECOND DISTANCE TO
C SOME LARGE
VALUE, SUCH AS 500000.0
C
DIPPROPS(1) = 1.296019
DIPPROPS(2) = 22.4717
DIPPROPS(3) = -7.22234
DIPPROPS(4) = 1465.76
DIPPROPS(5) = -838.57
C
C DIPPROPS ARE
THE DIPOLE PROPERTIES ALONG THE AXIS FROM AB
C INITIO
EVALUATION.
C DIPPROPS(1) = HCN DIPOLE
C
DIPPROPS(2) = HCN AXIAL DIPOLE POLAROZABILITY
C
AND SO ON
C
C FOLLOWING
ARE THE PARAMETERS FOR EQN. 6 OF THE PAPER BY
C C. E.
DYKSTRA
C
CHI15 =
2.0505478D 0
DEE15 =
0.0697675D 0
RRR15 =
1.341346D 0
DAA = 0.00056615D 0
DBB = 4.39881525D 0
CHIA = 0.5414377D 0
CHIB = 3.1438926D 0
SS1 = 0.737249D 0
SS2 = 1.004431D 0
SS3 = 1.355880D 0
SS4 = 2.261850D 0
MSTART =
0
C
10 PRINT*, ' ENTER THE PROTON-NITROGEN DISTANCES IN ANGSTROMS'
READ(5,*)
R1, R2
IF(MSTART
.EQ. 0) WRITE(6,940)
MSTART =
MSTART + 1
R12 = R1
+ R2
RR1 = (R1
+ 0.5950329) / 0.529177249D 0
RR2 = (R2
+ 0.5950329) / 0.529177249D 0
BIGR =
RR1 + RR2
BIGR3 =
BIGR ** 3
RR1RR1 =
RR1 * RR1
RR2RR2 =
RR2 * RR2
TWOMU =
2.0D 0 * DIPPROPS(1) / BIGR3
FIELD1 =
1.0D 0 / RR1RR1 - TWOMU
FIELD2 =
1.0D 0 / RR2RR2 - TWOMU
EELEC =
0.0D 0
FELEC =
0.0D 0
FACTOR =
1.0D 0
DO 100
N=1,5
FELEC =
FELEC - DIPPROPS(N) * FIELD1 ** N / FACTOR
FELEC =
FELEC - DIPPROPS(N) * FIELD2 ** N / FACTOR
EELEC =
EELEC - DIPPROPS(N) / (FACTOR * (RR1RR1 ** N))
EELEC =
EELEC - DIPPROPS(N) / (FACTOR * (RR2RR2 ** N))
FACTOR =
FACTOR * (N + 1.0D 0)
100 CONTINUE
ELEC3 =
FELEC - EELEC
ELECA =
2.0D 0 * DIPPROPS(2)**2 / (BIGR3 * RR1RR1 * RR2RR2)
ELECA = ELECA * (1.0D 0 - TWOMU * (RR1RR1 +
RR2RR2) )
ELEC3 =
ELEC3 + ELECA
EELEC =
FELEC + ELECA
WRITE(6,900) ELEC3, EELEC
C
C NEXT,
EVALUATE THE "PENETRATION" CONTRIBUTION (EQN. 6)
C
EPENET =
(1.0 - DEXP(-CHI15*(R1-RRR15)) ) ** 2
1
+(1.0 - DEXP(-CHI15*(R2-RRR15)) ) ** 2
EPENET =
DEE15 * (EPENET - 2.0D 0)
EELEC =
EELEC + EPENET
WRITE(6,910) EPENET, EELEC
C
C NEXT,
EVALUATE REFINEMENT VIA EQN. 7
C
EREFIN =
DAA * (R1 * R1 * DEXP(-CHIA * R1)
1
+R2 * R2 * DEXP(-CHIA * R2) )
2 +
DBB * ( (R1-SS1) * (R1-SS2) * (R1-SS3) * (R1-SS4)
3
* DEXP(-CHIB * R1)
4
+(R2-SS1) * (R2-SS2) * (R2-SS3) * (R2-SS4)
5
* DEXP(-CHIB * R2) )
EELEC =
EELEC + EREFIN
WRITE(6,920) EREFIN, EELEC
C
C
NEXT, ADD THE ISOLATED HCN-NCH INTERACTION ENERGY
C
(BASED ON A FIT OF AB INITIO VALUES)
C
EDIMER
= 0.390041956E-01 / (R12 ** 2)
1
-0.742043831 /
(R12 ** 6)
2 +
27.9499949 /
(R12 ** 8)
EELEC =
EELEC + EDIMER
WRITE(6,950) EDIMER, EELEC
C
C
NEXT, INCLUDE THE REMAINDER FITTING (NON-ADDITIVE PART)
C
RR1 = R1
+ 0.5950329D 0
RR2 = R2
+ 0.5950329D 0
R12 = R1
+ R2
DIFF2 =
(R1 - R2) ** 2
BIGR =
RR1 + RR2
BIGR3 =
BIGR ** 3
BIGR4 =
BIGR3 * BIGR
BIGR5 =
BIGR4 * BIGR
BIGR6 =
BIGR5 * BIGR
BIGR7 =
BIGR6 * BIGR
BIGR8 =
BIGR7 * BIGR
RRTT2 =
RR1 * RR1 * RR2 * RR2
EEEEE =
0.103355791 * EXP(-0.5*1.06516963*R12 *
1
(0.5*R12 - 1.163920999) )
EREM
= -46.8899545 /
RRTT2
1 + 163.770805
/ (RRTT2 * RR1 * RR2)
2 - 6.68441635
/ BIGR3
3 + 846.127015
/ BIGR4
4 - 448.231397
/ BIGR5
5 - 9806.99352
/ BIGR6
6 - 432.520703
* DIFF2 / BIGR4
7 + 7404.60412
* DIFF2 / BIGR5
8 - 42743.9909
* DIFF2 / BIGR6
9 + 105878.050
* DIFF2 / BIGR7
A - 114423.895
* DIFF2 / BIGR8
EREM =
EREM + EEEEE
EELEC =
EELEC + EREM
WRITE(6,930) EREM, EELEC, EELEC-EDIMER
GO TO 10
C
900 FORMAT(' 3-BODY AND NET ELECT. ENERGY IS ',2F18.9)
910 FORMAT(' THE PENETRATION ENERGY (A.U.) IS',2F18.9)
920 FORMAT(' THE REFINEMENT PENETRATION E IS ',2F18.9)
930 FORMAT(' REMAINING NON-ADDITIVE ENERGY IS',2F18.9/
1
2X,'('F14.9,' W/O
HCN-NCH)'/)
940 FORMAT(' INTERACTION ENERGY CONTRIBUTIONS IN A.U.'
1 '
ARE SUMMED IN THE LAST COLUMN')
950 FORMAT(' THE ISOLATED HCN-NCH ENERGY IS ',2F18.9)
END