! #+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+# ! #| Compute the B matrix |# ! #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#+#+#+#+# #include "../include/dprec.fh" subroutine wdc_bmat ( qcart, bmat ) use schlegel, only: ncoord, natm, ibond, iangle, idihed, nbond, nangle & , ndihed use evb_math, only: uvec, orthovec, vrot, dist implicit none _REAL_ , intent( in) :: qcart(3,natm) _REAL_ , intent(out) :: bmat(ncoord,natm*3) ! .......................................................................... integer :: i, j, k, l, m, n, jj, kk, ll, mm, nn _REAL_ :: s_tj(3), s_tk(3), s_tl(3), s_tm(3) bmat(:,:) = 0.0d0 nn = 0 ! +---------------------------------------------------------------------------+ ! | Compute the B matrix elements for a bond stretch | ! :...........................................................................: ! | Wilson, Decius, and Cross, "Molecular Vibrations" p. 55-56 | ! | | ! | j --> k r_jk = r_k - r_j | ! | | ! | s_tj = e_kj = - e_jk s_tk = e_jk = - e_kj | ! +---------------------------------------------------------------------------+ do n = 1, nbond j = ibond(n,1) k = ibond(n,2) s_tj(:) = uvec(qcart,k,j) jj = ( j - 1 ) * 3 kk = ( k - 1 ) * 3 nn = nn + 1 do i = 1, 3 bmat(nn,jj+i) = - s_tj(i) bmat(nn,kk+i) = s_tj(i) enddo enddo ! +---------------------------------------------------------------------------+ ! | Compute the B matrix elements for an angle bend | ! :...........................................................................: ! | Wilson, Decius, and Cross, "Molecular Vibrations" p. 56-58 | ! | | ! | l s_tj = ( cos(phi) e_lj - e_lk ) | ! | / \ / ( r_lj sin(phi) ) | ! | j k | ! | s_tk = ( cos(phi) e_lk - e_lj ) | ! | / ( r_lk sin(phi) ) | ! | | ! | s_tl = - s_tj - s_tk | ! | | ! | cos(phi) = e_lj * e_lk sin(phi) = SQRT( 1 - cos^2(phi) ) | ! +---------------------------------------------------------------------------+ do n = 1, nangle j = iangle(n,1) k = iangle(n,2) l = iangle(n,3) jj = ( j - 1 ) * 3 kk = ( k - 1 ) * 3 ll = ( l - 1 ) * 3 nn = nn + 1 s_tj(:) = - orthovec(qcart,k,l,j) / dist(qcart,j,l) s_tk(:) = - orthovec(qcart,j,l,k) / dist(qcart,k,l) s_tl(:) = - s_tj(:) - s_tk(:) do i = 1, 3 bmat(nn,jj+i) = s_tj(i) bmat(nn,kk+i) = s_tk(i) bmat(nn,ll+i) = s_tl(i) enddo enddo ! +---------------------------------------------------------------------------+ ! | Compute the B matrix elements for a proper dihedral | ! :...........................................................................: ! | Wilson, Decius, and Cross, "Molecular Vibrations" p. 60-61 | ! | | ! | s_tj = - e_jk X e_kl / r_jk / sin^2(phi2) | ! | j (phi2) | ! | \ s_tk = ( r_kl - r_jk cos(phi2) ) | ! | k -- l * e_jk X e_kl / r_kl / r_jk | ! | \ / sin^2(phi2) + cos(phi3) e_ml x e_lk | ! | (phi3) m / r_kl / sin^2(phi3) | ! | | ! | s_tl = [(jm)(kl)]_s_tk | ! | | ! | s_tm = [(jm)(kl)]_s_tj | ! | | ! | cos(tau) = (e_jk X e_kl) * (e_kl X e_lm) | ! | / sin(phi2) / sin(phi3) | ! | | ! | s_tl & s_tm are obtained via permutation of j with m and | ! | k with l and correponding phi2 and phi3 angles | ! +---------------------------------------------------------------------------+ do n = 1, ndihed j = idihed(n,1) k = idihed(n,2) l = idihed(n,3) m = idihed(n,4) jj = ( j - 1 ) * 3 kk = ( k - 1 ) * 3 ll = ( l - 1 ) * 3 mm = ( m - 1 ) * 3 nn = nn + 1 s_tj(:) = vrot(qcart,j,k,l) s_tk(:) = - vrot(qcart,j,k,l) * ( 1.0d0 - dot_product( uvec(qcart,j,k) & , uvec(qcart,l,k) ) * dist(qcart,j,k) / dist(qcart,k,l) ) & - vrot(qcart,m,l,k) * dot_product( uvec(qcart,k,l) & , uvec(qcart,m,l) ) * dist(qcart,l,m) / dist(qcart,k,l) s_tl(:) = - vrot(qcart,m,l,k) * ( 1.0d0 - dot_product( uvec(qcart,k,l) & , uvec(qcart,m,l) ) * dist(qcart,l,m) / dist(qcart,k,l) ) & - vrot(qcart,j,k,l) * dot_product( uvec(qcart,j,k) & , uvec(qcart,l,k) ) * dist(qcart,j,k) / dist(qcart,k,l) s_tm(:) = vrot(qcart,m,l,k) do i = 1, 3 bmat(nn,jj+i) = s_tj(i) bmat(nn,kk+i) = s_tk(i) bmat(nn,ll+i) = s_tl(i) bmat(nn,mm+i) = s_tm(i) enddo enddo ! x+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++x ! x+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++x end subroutine wdc_bmat