1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ [Enter a one-line description of subroutine difbon here]
subroutine difbon(dd,df,ddf,xij,b, &
i3p,j3p,nbel)
implicit double precision (a-h,o-z)
! ----- calculate second derivatives of bond-like terms -----
! dd ... second derivative matrix
! df ... first der. of v with respect to bond length
! ddf ... corresponding second der.
! xij ... bond vector x(i)-x(j) etc
! b ... bond length
! i3p ... atom coordinate index for atom i (3*i-3)
! j3p ... atom coordinate for atom j
! coded by d.a.case
dimension dr(6),ddr(6,6), nbel(*), istore(6,6)
dimension dd(2),xij(2)
ibel(m3) = nbel(m3/3+1)
i3 = ibel(min(i3p,j3p))
j3 = ibel(max(i3p,j3p))
! ----- first set up array dr, which holds first derivative of bond
! length with respect to cartesian coordinates -----
do 10 k = 1,3
dr(k) = xij(k)/b
dr(k+3) = -dr(k)
10 continue
! ----- now set up array ddr, which contains the second derivatives
! of bond length with respect to cartesians -----
a = (1.0e0-dr(1)*dr(1))/b
ddr(1,1) = a
ddr(4,4) = a
ddr(1,4) = -a
a = (1.0e0-dr(2)*dr(2))/b
ddr(2,2) = a
ddr(5,5) = a
ddr(2,5) = -a
a = (1.0e0-dr(3)*dr(3))/b
ddr(3,3) = a
ddr(6,6) = a
ddr(3,6) = -a
a = -dr(1)*dr(2)/b
ddr(1,2) = a
ddr(4,5) = a
ddr(1,5) = -a
ddr(2,4) = -a
a = -dr(1)*dr(3)/b
ddr(1,3) = a
ddr(4,6) = a
ddr(1,6) = -a
ddr(3,4) = -a
a = -dr(2)*dr(3)/b
ddr(2,3) = a
ddr(5,6) = a
ddr(2,6) = -a
ddr(3,5) = -a
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!******************************************************************************************
!* INSERTION DES INSTRUCTIONS POUR L'ECRITURE DANS UN FICHIER deriv.txt DE LA MATRICE ddr *
!******************************************************************************************
OPEN(unit=10, file='deriv.txt', iostat=ios)
IF (IOS /= 0) THEN
WRITE(6, *) 'Could not create deriv.txt'
! gerer l'erreur -> fin du programme (violent) ou retour d'un code
! d'erreur a la fonction principale
RETURN ios ! ou STOP
ENDIF
DO column = 1, 6
DO line = 1, 6
! ici, on peut faire le calcul de values(i, j)
! on ecrit la valeur dans le fichier ouvert precedemment
WRITE(10, fmt='(1pe12.5)') ddr(line,column),'|'
ENDDO
write(unit=1, fmt=*) 'newline'
ENDDO
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
