TABLE OF CONTENTS
ABINIT/cont33so [ Functions ]
NAME
cont33so
FUNCTION
Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor gxa2 using metric tensor gmet and antisymmetric tensor amet to produce rank 2 real tensor.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
gxa1(2,10)=rank 3 complex symmetric tensor gxa2(2,10)=rank 3 complex symmetric tensor gmet(3,3)=usual metric tensor (symmetric, real) amet(2,3,3)=antisymmetric complex tensor used for spin-orbit
OUTPUT
rank2(6)=rank 2 real tensor (pseudo-symmetric storage)
NOTES
This contraction is used for spin-orbit correction in non-local contribution to stresses. Symmetric gxa1, gxa2 are stored as 111 221 331 321 311 211 222 332 322 333; gmet(3,3) is symmetric but stored fully (9 elements); amet(3,3) is antisymmetric but stored fully (9 elements); Output rank2 is not symmetric but since $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$ it is stored as 11 22 33 32 31 21. Want 2*Re[contraction]. rank2(a,b)=2 Re[15 r_{3A}(a,i,j) r_{3G}(b,j,i)-3 r_{3A}(a,b,i) r_1(i)] where: r_1(i) & = & gxa2(i,j,k) gmet(j,k) \nonumber r_{3A}(i,j,k) & = & conjg(gxa1(p,i,j)) amet(p,k) \nonumber r_{3G}(i,j,k) & = & gxa2(p,i,j) gmet(p,k) \end{eqnarray} }}
PARENTS
nonlop_pl
CHILDREN
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 64 subroutine cont33so(gxa1,gxa2,gmet,amet,rank2) 65 66 use defs_basis 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'cont33so' 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !arrays 78 real(dp),intent(in) :: amet(2,3,3),gmet(3,3),gxa1(2,10),gxa2(2,10) 79 real(dp),intent(out) :: rank2(6) 80 81 !Local variables------------------------------- 82 !scalars 83 integer,parameter :: im=2,re=1 84 integer :: ii 85 !arrays 86 real(dp) :: r1(2,3),r3A(2,18),r3G(2,18),s31(6),s33(6) 87 88 ! ************************************************************************* 89 90 !Compute r1(i)=gxa2(i,j,k)*gmet(j,k) 91 92 do ii=1,2 93 r1(ii,1)=gmet(1,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,2)+& 94 & gmet(3,3)*gxa2(ii,3)+2.d0*(& 95 & gmet(3,2)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,5)+& 96 & gmet(2,1)*gxa2(ii,6)) 97 r1(ii,2)=gmet(1,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,7)+& 98 & gmet(3,3)*gxa2(ii,8)+2.d0*(& 99 & gmet(3,2)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,4)+& 100 & gmet(2,1)*gxa2(ii,2)) 101 r1(ii,3)=gmet(1,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,9)+& 102 & gmet(3,3)*gxa2(ii,10)+2.d0*(& 103 & gmet(3,2)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,3)+& 104 & gmet(2,1)*gxa2(ii,4)) 105 end do 106 107 !Compute r3G(i,j,k)=gxa2(p,i,j)*gmet(p,k) 108 !Write out components for 18 distinct terms, Re and Im 109 !(r3G is symmetric in first two indices, not in all permutations) 110 !Store as 111 221 331 321 311 211 111 !112 222 332 322 312 212 112 !113 223 333 323 313 213 113 do ii=1,2 114 r3G(ii, 1)=gmet(1,1)*gxa2(ii,1)+gmet(2,1)*gxa2(ii,6)+gmet(3,1)*gxa2(ii,5) 115 r3G(ii, 2)=gmet(1,1)*gxa2(ii,2)+gmet(2,1)*gxa2(ii,7)+gmet(3,1)*gxa2(ii,9) 116 r3G(ii, 3)=gmet(1,1)*gxa2(ii,3)+gmet(2,1)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,10) 117 r3G(ii, 4)=gmet(1,1)*gxa2(ii,4)+gmet(2,1)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,8) 118 r3G(ii, 5)=gmet(1,1)*gxa2(ii,5)+gmet(2,1)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,3) 119 r3G(ii, 6)=gmet(1,1)*gxa2(ii,6)+gmet(2,1)*gxa2(ii,2)+gmet(3,1)*gxa2(ii,4) 120 r3G(ii, 7)=gmet(2,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,5) 121 r3G(ii, 8)=gmet(2,1)*gxa2(ii,2)+gmet(2,2)*gxa2(ii,7)+gmet(3,2)*gxa2(ii,9) 122 r3G(ii, 9)=gmet(2,1)*gxa2(ii,3)+gmet(2,2)*gxa2(ii,8)+gmet(3,2)*gxa2(ii,10) 123 r3G(ii,10)=gmet(2,1)*gxa2(ii,4)+gmet(2,2)*gxa2(ii,9)+gmet(3,2)*gxa2(ii,8) 124 r3G(ii,11)=gmet(2,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,3) 125 r3G(ii,12)=gmet(2,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,4) 126 r3G(ii,13)=gmet(3,1)*gxa2(ii,1)+gmet(3,2)*gxa2(ii,6)+gmet(3,3)*gxa2(ii,5) 127 r3G(ii,14)=gmet(3,1)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,7)+gmet(3,3)*gxa2(ii,9) 128 r3G(ii,15)=gmet(3,1)*gxa2(ii,3)+gmet(3,2)*gxa2(ii,8)+gmet(3,3)*gxa2(ii,10) 129 r3G(ii,16)=gmet(3,1)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,9)+gmet(3,3)*gxa2(ii,8) 130 r3G(ii,17)=gmet(3,1)*gxa2(ii,5)+gmet(3,2)*gxa2(ii,4)+gmet(3,3)*gxa2(ii,3) 131 r3G(ii,18)=gmet(3,1)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,2)+gmet(3,3)*gxa2(ii,4) 132 end do 133 134 !Compute r3A(i,j,k)=conjg(gxa1(p,i,j))*amet(p,k) 135 !Write out components for 18 distinct terms, Re and Im 136 !(r3A is symmetric in first two indices, not in all permutations) 137 !Store as 111 221 331 321 311 211 138 !112 222 332 322 312 212 139 !113 223 333 323 313 213 140 !Note that, since amet is antisymmetric, amet(i,i)=0 141 142 r3A(re, 1)=amet(re,2,1)*gxa1(re,6 )+amet(im,2,1)*gxa1(im,6 )+& 143 & amet(re,3,1)*gxa1(re,5 )+amet(im,3,1)*gxa1(im,5 ) 144 r3A(re, 2)=amet(re,2,1)*gxa1(re,7 )+amet(im,2,1)*gxa1(im,7 )+& 145 & amet(re,3,1)*gxa1(re,9 )+amet(im,3,1)*gxa1(im,9 ) 146 r3A(re, 3)=amet(re,2,1)*gxa1(re,8 )+amet(im,2,1)*gxa1(im,8 )+& 147 & amet(re,3,1)*gxa1(re,10)+amet(im,3,1)*gxa1(im,10) 148 r3A(re, 4)=amet(re,2,1)*gxa1(re,9 )+amet(im,2,1)*gxa1(im,9 )+& 149 & amet(re,3,1)*gxa1(re,8 )+amet(im,3,1)*gxa1(im,8 ) 150 r3A(re, 5)=amet(re,2,1)*gxa1(re,4 )+amet(im,2,1)*gxa1(im,4 )+& 151 & amet(re,3,1)*gxa1(re,3 )+amet(im,3,1)*gxa1(im,3 ) 152 r3A(re, 6)=amet(re,2,1)*gxa1(re,2 )+amet(im,2,1)*gxa1(im,2 )+& 153 & amet(re,3,1)*gxa1(re,4 )+amet(im,3,1)*gxa1(im,4 ) 154 r3A(re, 7)=amet(re,1,2)*gxa1(re,1 )+amet(im,1,2)*gxa1(im,1 )+& 155 & amet(re,3,2)*gxa1(re,5 )+amet(im,3,2)*gxa1(im,5 ) 156 r3A(re, 8)=amet(re,1,2)*gxa1(re,2 )+amet(im,1,2)*gxa1(im,2 )+& 157 & amet(re,3,2)*gxa1(re,9 )+amet(im,3,2)*gxa1(im,9 ) 158 r3A(re, 9)=amet(re,1,2)*gxa1(re,3 )+amet(im,1,2)*gxa1(im,3 )+& 159 & amet(re,3,2)*gxa1(re,10)+amet(im,3,2)*gxa1(im,10) 160 r3A(re,10)=amet(re,1,2)*gxa1(re,4 )+amet(im,1,2)*gxa1(im,4 )+& 161 & amet(re,3,2)*gxa1(re,8 )+amet(im,3,2)*gxa1(im,8 ) 162 r3A(re,11)=amet(re,1,2)*gxa1(re,5 )+amet(im,1,2)*gxa1(im,5 )+& 163 & amet(re,3,2)*gxa1(re,3 )+amet(im,3,2)*gxa1(im,3 ) 164 r3A(re,12)=amet(re,1,2)*gxa1(re,6 )+amet(im,1,2)*gxa1(im,6 )+& 165 & amet(re,3,2)*gxa1(re,4 )+amet(im,3,2)*gxa1(im,4 ) 166 r3A(re,13)=amet(re,1,3)*gxa1(re,1 )+amet(im,1,3)*gxa1(im,1 )+& 167 & amet(re,2,3)*gxa1(re,6 )+amet(im,2,3)*gxa1(im,6 ) 168 r3A(re,14)=amet(re,1,3)*gxa1(re,2 )+amet(im,1,3)*gxa1(im,2 )+& 169 & amet(re,2,3)*gxa1(re,7 )+amet(im,2,3)*gxa1(im,7 ) 170 r3A(re,15)=amet(re,1,3)*gxa1(re,3 )+amet(im,1,3)*gxa1(im,3 )+& 171 & amet(re,2,3)*gxa1(re,8 )+amet(im,2,3)*gxa1(im,8 ) 172 r3A(re,16)=amet(re,1,3)*gxa1(re,4 )+amet(im,1,3)*gxa1(im,4 )+& 173 & amet(re,2,3)*gxa1(re,9 )+amet(im,2,3)*gxa1(im,9 ) 174 r3A(re,17)=amet(re,1,3)*gxa1(re,5 )+amet(im,1,3)*gxa1(im,5 )+& 175 & amet(re,2,3)*gxa1(re,4 )+amet(im,2,3)*gxa1(im,4 ) 176 r3A(re,18)=amet(re,1,3)*gxa1(re,6 )+amet(im,1,3)*gxa1(im,6 )+& 177 & amet(re,2,3)*gxa1(re,2 )+amet(im,2,3)*gxa1(im,2 ) 178 179 r3A(im, 1)=amet(im,2,1)*gxa1(re,6 )-amet(re,2,1)*gxa1(im,6 )+& 180 & amet(im,3,1)*gxa1(re,5 )-amet(re,3,1)*gxa1(im,5 ) 181 r3A(im, 2)=amet(im,2,1)*gxa1(re,7 )-amet(re,2,1)*gxa1(im,7 )+& 182 & amet(im,3,1)*gxa1(re,9 )-amet(re,3,1)*gxa1(im,9 ) 183 r3A(im, 3)=amet(im,2,1)*gxa1(re,8 )-amet(re,2,1)*gxa1(im,8 )+& 184 & amet(im,3,1)*gxa1(re,10)-amet(re,3,1)*gxa1(im,10) 185 r3A(im, 4)=amet(im,2,1)*gxa1(re,9 )-amet(re,2,1)*gxa1(im,9 )+& 186 & amet(im,3,1)*gxa1(re,8 )-amet(re,3,1)*gxa1(im,8 ) 187 r3A(im, 5)=amet(im,2,1)*gxa1(re,4 )-amet(re,2,1)*gxa1(im,4 )+& 188 & amet(im,3,1)*gxa1(re,3 )-amet(re,3,1)*gxa1(im,3 ) 189 r3A(im, 6)=amet(im,2,1)*gxa1(re,2 )-amet(re,2,1)*gxa1(im,2 )+& 190 & amet(im,3,1)*gxa1(re,4 )-amet(re,3,1)*gxa1(im,4 ) 191 r3A(im, 7)=amet(im,1,2)*gxa1(re,1 )-amet(re,1,2)*gxa1(im,1 )+& 192 & amet(im,3,2)*gxa1(re,5 )-amet(re,3,2)*gxa1(im,5 ) 193 r3A(im, 8)=amet(im,1,2)*gxa1(re,2 )-amet(re,1,2)*gxa1(im,2 )+& 194 & amet(im,3,2)*gxa1(re,9 )-amet(re,3,2)*gxa1(im,9 ) 195 r3A(im, 9)=amet(im,1,2)*gxa1(re,3 )-amet(re,1,2)*gxa1(im,3 )+& 196 & amet(im,3,2)*gxa1(re,10)-amet(re,3,2)*gxa1(im,10) 197 r3A(im,10)=amet(im,1,2)*gxa1(re,4 )-amet(re,1,2)*gxa1(im,4 )+& 198 & amet(im,3,2)*gxa1(re,8 )-amet(re,3,2)*gxa1(im,8 ) 199 r3A(im,11)=amet(im,1,2)*gxa1(re,5 )-amet(re,1,2)*gxa1(im,5 )+& 200 & amet(im,3,2)*gxa1(re,3 )-amet(re,3,2)*gxa1(im,3 ) 201 r3A(im,12)=amet(im,1,2)*gxa1(re,6 )-amet(re,1,2)*gxa1(im,6 )+& 202 & amet(im,3,2)*gxa1(re,4 )-amet(re,3,2)*gxa1(im,4 ) 203 r3A(im,13)=amet(im,1,3)*gxa1(re,1 )-amet(re,1,3)*gxa1(im,1 )+& 204 & amet(im,2,3)*gxa1(re,6 )-amet(re,2,3)*gxa1(im,6 ) 205 r3A(im,14)=amet(im,1,3)*gxa1(re,2 )-amet(re,1,3)*gxa1(im,2 )+& 206 & amet(im,2,3)*gxa1(re,7 )-amet(re,2,3)*gxa1(im,7 ) 207 r3A(im,15)=amet(im,1,3)*gxa1(re,3 )-amet(re,1,3)*gxa1(im,3 )+& 208 & amet(im,2,3)*gxa1(re,8 )-amet(re,2,3)*gxa1(im,8 ) 209 r3A(im,16)=amet(im,1,3)*gxa1(re,4 )-amet(re,1,3)*gxa1(im,4 )+& 210 & amet(im,2,3)*gxa1(re,9 )-amet(re,2,3)*gxa1(im,9 ) 211 r3A(im,17)=amet(im,1,3)*gxa1(re,5 )-amet(re,1,3)*gxa1(im,5 )+& 212 & amet(im,2,3)*gxa1(re,4 )-amet(re,2,3)*gxa1(im,4 ) 213 r3A(im,18)=amet(im,1,3)*gxa1(re,6 )-amet(re,1,3)*gxa1(im,6 )+& 214 & amet(im,2,3)*gxa1(re,2 )-amet(re,2,3)*gxa1(im,2 ) 215 216 !Compute s33(a,b)=2*Re[r3A(a,i,j)*r3G(b,j,i)] 217 218 s33(1)=2.d0*(r3A(re, 1)*r3G(re, 1)-r3A(im, 1)*r3G(im, 1)+& 219 & r3A(re,12)*r3G(re,12)-r3A(im,12)*r3G(im,12)+& 220 & r3A(re,17)*r3G(re,17)-r3A(im,17)*r3G(im,17)+& 221 & r3A(re,11)*r3G(re,18)-r3A(im,11)*r3G(im,18)+& 222 & r3A(re,18)*r3G(re,11)-r3A(im,18)*r3G(im,11)+& 223 & r3A(re, 5)*r3G(re,13)-r3A(im, 5)*r3G(im,13)+& 224 & r3A(re,13)*r3G(re, 5)-r3A(im,13)*r3G(im, 5)+& 225 & r3A(re, 6)*r3G(re, 7)-r3A(im, 6)*r3G(im, 7)+& 226 & r3A(re, 7)*r3G(re, 6)-r3A(im, 7)*r3G(im, 6)) 227 s33(2)=2.d0*(r3A(re, 6)*r3G(re, 6)-r3A(im, 6)*r3G(im, 6)+& 228 & r3A(re, 8)*r3G(re, 8)-r3A(im, 8)*r3G(im, 8)+& 229 & r3A(re,16)*r3G(re,16)-r3A(im,16)*r3G(im,16)+& 230 & r3A(re,10)*r3G(re,14)-r3A(im,10)*r3G(im,14)+& 231 & r3A(re,14)*r3G(re,10)-r3A(im,14)*r3G(im,10)+& 232 & r3A(re, 4)*r3G(re,18)-r3A(im, 4)*r3G(im,18)+& 233 & r3A(re,18)*r3G(re, 4)-r3A(im,18)*r3G(im, 4)+& 234 & r3A(re, 2)*r3G(re,12)-r3A(im, 2)*r3G(im,12)+& 235 & r3A(re,12)*r3G(re, 2)-r3A(im,12)*r3G(im, 2)) 236 s33(3)=2.d0*(r3A(re, 5)*r3G(re, 5)-r3A(im, 5)*r3G(im, 5)+& 237 & r3A(re,10)*r3G(re,10)-r3A(im,10)*r3G(im,10)+& 238 & r3A(re,15)*r3G(re,15)-r3A(im,15)*r3G(im,15)+& 239 & r3A(re, 9)*r3G(re,16)-r3A(im, 9)*r3G(im,16)+& 240 & r3A(re,16)*r3G(re, 9)-r3A(im,16)*r3G(im, 9)+& 241 & r3A(re, 3)*r3G(re,17)-r3A(im, 3)*r3G(im,17)+& 242 & r3A(re,17)*r3G(re, 3)-r3A(im,17)*r3G(im, 3)+& 243 & r3A(re, 4)*r3G(re,11)-r3A(im, 4)*r3G(im,11)+& 244 & r3A(re,11)*r3G(re, 4)-r3A(im,11)*r3G(im, 4)) 245 s33(4)=2.d0*(r3A(re, 5)*r3G(re, 6)-r3A(im, 5)*r3G(im, 6)+& 246 & r3A(re,10)*r3G(re, 8)-r3A(im,10)*r3G(im, 8)+& 247 & r3A(re,15)*r3G(re,16)-r3A(im,15)*r3G(im,16)+& 248 & r3A(re, 9)*r3G(re,14)-r3A(im, 9)*r3G(im,14)+& 249 & r3A(re,16)*r3G(re,10)-r3A(im,16)*r3G(im,10)+& 250 & r3A(re, 3)*r3G(re,18)-r3A(im, 3)*r3G(im,18)+& 251 & r3A(re,17)*r3G(re, 4)-r3A(im,17)*r3G(im, 4)+& 252 & r3A(re, 4)*r3G(re,12)-r3A(im, 4)*r3G(im,12)+& 253 & r3A(re,11)*r3G(re, 2)-r3A(im,11)*r3G(im, 2)) 254 s33(5)=2.d0*(r3A(re, 5)*r3G(re, 1)-r3A(im, 5)*r3G(im, 1)+& 255 & r3A(re,10)*r3G(re,12)-r3A(im,10)*r3G(im,12)+& 256 & r3A(re,15)*r3G(re,17)-r3A(im,15)*r3G(im,17)+& 257 & r3A(re, 9)*r3G(re,18)-r3A(im, 9)*r3G(im,18)+& 258 & r3A(re,16)*r3G(re,11)-r3A(im,16)*r3G(im,11)+& 259 & r3A(re, 3)*r3G(re,13)-r3A(im, 3)*r3G(im,13)+& 260 & r3A(re,17)*r3G(re, 5)-r3A(im,17)*r3G(im, 5)+& 261 & r3A(re, 4)*r3G(re, 7)-r3A(im, 4)*r3G(im, 7)+& 262 & r3A(re,11)*r3G(re, 6)-r3A(im,11)*r3G(im, 6)) 263 s33(6)=2.d0*(r3A(re, 6)*r3G(re, 1)-r3A(im, 6)*r3G(im, 1)+& 264 & r3A(re, 8)*r3G(re,12)-r3A(im, 8)*r3G(im,12)+& 265 & r3A(re,16)*r3G(re,17)-r3A(im,16)*r3G(im,17)+& 266 & r3A(re,10)*r3G(re,18)-r3A(im,10)*r3G(im,18)+& 267 & r3A(re,14)*r3G(re,11)-r3A(im,14)*r3G(im,11)+& 268 & r3A(re, 4)*r3G(re,13)-r3A(im, 4)*r3G(im,13)+& 269 & r3A(re,18)*r3G(re, 5)-r3A(im,18)*r3G(im, 5)+& 270 & r3A(re, 2)*r3G(re, 7)-r3A(im, 2)*r3G(im, 7)+& 271 & r3A(re,12)*r3G(re, 6)-r3A(im,12)*r3G(im, 6)) 272 273 !Compute s31(a,b)=2*Re[r3A(a,b,i)*r1(i)] 274 275 s31(1)=2.d0*(r1(re,1)*r3A(re, 1)-r1(im,1)*r3A(im, 1)+& 276 & r1(re,2)*r3A(re, 7)-r1(im,2)*r3A(im, 7)+& 277 & r1(re,3)*r3A(re,13)-r1(im,3)*r3A(im,13)) 278 s31(2)=2.d0*(r1(re,1)*r3A(re, 2)-r1(im,1)*r3A(im, 2)+& 279 & r1(re,2)*r3A(re, 8)-r1(im,2)*r3A(im, 8)+& 280 & r1(re,3)*r3A(re,14)-r1(im,3)*r3A(im,14)) 281 s31(3)=2.d0*(r1(re,1)*r3A(re, 3)-r1(im,1)*r3A(im, 3)+& 282 & r1(re,2)*r3A(re, 9)-r1(im,2)*r3A(im, 9)+& 283 & r1(re,3)*r3A(re,15)-r1(im,3)*r3A(im,15)) 284 s31(4)=2.d0*(r1(re,1)*r3A(re, 4)-r1(im,1)*r3A(im, 4)+& 285 & r1(re,2)*r3A(re,10)-r1(im,2)*r3A(im,10)+& 286 & r1(re,3)*r3A(re,16)-r1(im,3)*r3A(im,16)) 287 s31(5)=2.d0*(r1(re,1)*r3A(re, 5)-r1(im,1)*r3A(im, 5)+& 288 & r1(re,2)*r3A(re,11)-r1(im,2)*r3A(im,11)+& 289 & r1(re,3)*r3A(re,17)-r1(im,3)*r3A(im,17)) 290 s31(6)=2.d0*(r1(re,1)*r3A(re, 6)-r1(im,1)*r3A(im, 6)+& 291 & r1(re,2)*r3A(re,12)-r1(im,2)*r3A(im,12)+& 292 & r1(re,3)*r3A(re,18)-r1(im,3)*r3A(im,18)) 293 294 !Finally, compute rank2(a,b)=-15*s33(a,b)+3*s31(a,b) 295 296 rank2(:)=15.d0*s33(:)-3.d0*s31(:) 297 298 end subroutine cont33so