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
93
|
--- a/DALTON/soppa/so_stoppw.F
+++ b/DALTON/soppa/so_stoppw.F
@@ -1,19 +1,22 @@
C
-C /* Deck so_stoppw */
- SUBROUTINE SO_STOPPW(STOPP,TRGOS,ISYMTR,IEXCI,EXENG,QVAL)
+ SUBROUTINE SO_STOPPW(RSTOPP,TRGOS,ISYMTR,IEXCI,EXENG,QVAL)
C
C This routine is part of the atomic integral directSOPPA program.
-C
+C The charge Z of the incoming ions is set to 1 here.
C Zhiwen Shi (Clark), Stephan P. A. Sauer, January 2016
-C
C PURPOSE: Calculate Stopping Power.
C
-#include "implicit.h"
-#include "cbiexc.h"
-#include "ccorb.h"
+ implicit none
+#include "cbiexc.h" ! LVEL, MXNEXI
+#include "ccorb.h" ! NSYM
+#include "pi.h" ! PI
C
- DIMENSION STOPP(3,LVEL),TRGOS(3),EXENG(NSYM,MXNEXI)
- REAL*8 QVAL,VELOC,QMAXV,QMINV
+ REAL*8, INTENT(INOUT) :: RSTOPP(3,LVEL,2)
+ REAL*8, INTENT(IN) :: TRGOS(3), EXENG(NSYM,MXNEXI), QVAL
+ INTEGER, INTENT(IN) :: ISYMTR, IEXCI
+ REAL*8 :: VELOC, QMAXV, QMINV
+ REAL*8, PARAMETER :: D4 = 4.0D0
+ INTEGER :: IVEL
C
C------------------
C Add to trace.
@@ -28,19 +31,52 @@ C
DO IVEL=1, LVEL
C
VELOC = VMIN+(IVEL-1)*VSTEP
-C
QMAXV = VELOC*2
QMINV = EXENG(ISYMTR,IEXCI)/VELOC
C
IF (QMINV .LE. QMAXV) THEN
C
- IF ((QVAL .GE. QMINV) .AND. (QVAL .LE. QMAXV)) THEN
+ IF ((QVAL .GE. QMINV) .AND. (QVAL .LE. QMAXV)) THEN
+C
+ RSTOPP(1,IVEL,1) = RSTOPP(1,IVEL,1) + TRGOS(1)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(2,IVEL,1) = RSTOPP(2,IVEL,1) + TRGOS(2)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(3,IVEL,1) = RSTOPP(3,IVEL,1) + TRGOS(3)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+
+ ENDIF
+C
+C The next line makes sure that the split integration will be applied
+C after the velocity is larger than half of Q value.
+C This velocity as half of chosen Q value corresponds to
+C the highest excitation energy for given basis set.
+C i.e. integration can be split after this velocity.
+C
+ IF (VELOC .GE. QINP/2.0d0) THEN
+C
+ IF ((QVAL .GE. QMINV) .AND. (QVAL .LE. QINP)
+ & .AND. (QVAL .LE. QMAXV)) THEN
+C
+ RSTOPP(1,IVEL,2) = RSTOPP(1,IVEL,2) + TRGOS(1)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(2,IVEL,2) = RSTOPP(2,IVEL,2) + TRGOS(2)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(3,IVEL,2) = RSTOPP(3,IVEL,2) + TRGOS(3)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+C
+ ENDIF
+C
+ ELSEIF ((QVAL .GE. QMINV) .AND. (QVAL .LE. QMAXV)) THEN
C
- STOPP(1,IVEL) = STOPP(1,IVEL) + TRGOS(1)*QSTEP
- STOPP(2,IVEL) = STOPP(2,IVEL) + TRGOS(2)*QSTEP
- STOPP(3,IVEL) = STOPP(3,IVEL) + TRGOS(3)*QSTEP
+ RSTOPP(1,IVEL,2) = RSTOPP(1,IVEL,2) + TRGOS(1)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(2,IVEL,2) = RSTOPP(2,IVEL,2) + TRGOS(2)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
+ RSTOPP(3,IVEL,2) = RSTOPP(3,IVEL,2) + TRGOS(3)*QSTEP*
+ & PI*D4/(QVAL*VELOC*VELOC)
C
- ENDIF
+ ENDIF
C
ENDIF
C
|
