Advertisement
Guest User

forces1

a guest
Nov 24th, 2017
86
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 6.51 KB | None | 0 0
  1. Module force_cal1
  2. use types
  3. use global_variables
  4. use math_subroutines
  5. use neigbours1
  6. use potpar1
  7. implicit none
  8. save
  9.  
  10.  
  11.  
  12. Contains
  13.  
  14. subroutine force_calc(v,xc,deriv)
  15. real(rk)::v,xc(3,n_at)
  16. real(rk)::epsmix,sigmix,sig2,sig6,sig12,vij,r1,fij(3)
  17. real(rk)::vcmix,dvcmix,dvccmix
  18. integer(ik)::i,j,k,ndev(3)
  19. logical::deriv
  20.  
  21.  
  22.  
  23.  
  24. if (pottyp.eq.0) then
  25. V=0.0_rk
  26. if (simtyp.eq.1) force=0.0_rk
  27. endif
  28. if (pottyp.eq.1) then
  29. rvec=0.0_rk
  30. v=0.0_rk
  31. force=0.0_rk
  32. P=0.0_rk
  33. der_boxlen=0.0_rk
  34. do i=1,n_at-1
  35. do j=i+1,n_at
  36. !r2=rsq(xc(:,i),xc(:,j),boxlen)
  37. call rsq_cal(xc(:,i),xc(:,j),boxlen,r2,rvec,ndev)
  38. v=r2**(-6)
  39. fij=2*(xc(:,i)-xc(:,j))/(r2**7)
  40. force(:,i)=force(:,i)+fij
  41. force(:,j)=force(:,j)-fij
  42. enddo
  43. enddo
  44. end if
  45. if (pottyp.eq.2) then
  46. call numderv2
  47.  
  48. rvec=0.0_rk
  49. v=0.0_rk
  50. force=0.0_rk
  51. P=0.0_rk
  52. der_boxlen=0.0_rk
  53. do i=1,n_at-1
  54. do j=i+1,n_at
  55. r2=rsq(xc(:,i),xc(:,j),boxlen)
  56. v=r2**(-1)
  57. fij=2*(xc(:,i)-xc(:,j))/(r2**2)
  58. force(:,i)=force(:,i)+fij
  59. force(:,j)=force(:,j)-fij
  60. !derivadas
  61. write(*,*) "derivadas nuemricas"
  62. ! call numder(valor,2)
  63. write(*,*) "derivadas analitica"
  64. write(*,*) force(:,1)
  65. write(*,*) force(:,2)
  66. stop
  67. enddo
  68. enddo
  69. !potencial igual a 1 y 3 con esfera rigida y esfera blanda
  70. end if
  71. !Escribir potenciales y fuerzas para esfera blanda y rígida
  72. if (pottyp.eq.3) then
  73. rvec=0.0_rk
  74. v=0.0_rk
  75. force=0.0_rk
  76. P=0.0_rk
  77. der_boxlen=0.0_rk
  78. do i=1,n_at-1
  79. do j=i+1,n_at
  80. call rsq_cal(xc(:,i),xc(:,j),boxlen,r2,rvec,ndev)
  81. if (r2.gt.rc2) cycle
  82. r1=dsqrt(r2)
  83. if (att(i)%atnum.ne.att(j)%atnum) then
  84. call mixpar(i,j,sigmix,epsmix)
  85. vcmix=max(lj(i)%vc,lj(j)%vc)
  86. dvcmix=max(lj(i)%dvc,lj(j)%dvc)
  87. dvccmix=max(lj(i)%dvcc,lj(j)%dvcc)
  88.  
  89. else
  90. epsmix=lj(i)%epslon
  91. sigmix=lj(i)%sigma
  92. vcmix=lj(i)%vc
  93. dvcmix=lj(i)%dvc
  94. dvccmix=lj(i)%dvcc
  95. endif
  96.  
  97. sig2=sigmix**2/r2
  98. sig6=sig2*sig2*sig2
  99. sig12=sig6*sig6
  100. v=v+4.0_rk*epsmix*(sig12-sig6-vcmix+dvccmix*(r1-rcut))
  101. if (deriv) then
  102. fij=24.0_rk*epsmix*(2*sig12-sig6-dvcmix*r1)*rvec/r2
  103. P=P+dot_product(fij,rvec)
  104. force(:,i)=force(:,i)+fij
  105. force(:,j)=force(:,j)-fij
  106. do k=1,3
  107. if (ndev(k).ne.0) then
  108. der_boxlen(k)=der_boxlen(k)+ndev(k)*fij(k)
  109. endif
  110. enddo
  111. endif
  112.  
  113. enddo
  114. enddo
  115.  
  116. endif
  117.  
  118. return
  119. end subroutine force_calc
  120.  
  121.  
  122.  
  123. !subroutine graph
  124. !integer(ik)::i
  125. !real(rk)::sig2,sig6,sig12,v,rvec(3),r2,fij(3),r
  126.  
  127.  
  128. !do i=0,1200
  129. ! r=0.5_rk+0.01_rk*i
  130. ! r2=r*r
  131. ! rvec=(/0.0_rk,0.0_rk,r/)
  132. ! v=0.0_rk
  133. ! sig2=sigma**2/r2
  134. ! sig6=sig2*sig2*sig2
  135. ! sig12=sig6*sig6
  136. ! if (r2.le.rc2) then
  137. ! v=v+sig12-sig6-vc+dvcc*(r-rcut)
  138. ! fij=(2.0_rk*sig12-sig6-dvc*r)*rvec/r2
  139. ! v=4.0_rk*epslon*v
  140. ! fij=24.0_rk*epslon*v
  141. ! endif
  142. ! write(1334,*) r,V,fij(3)
  143. !enddo
  144.  
  145.  
  146.  
  147.  
  148.  
  149. !end subroutine graph
  150.  
  151. subroutine mixpar(i,j,sigmix,epsmix)
  152. integer(ik),intent(in)::i,j
  153. real(rk)::sigmix,epsmix
  154.  
  155.  
  156. sigmix=0.5_rk*(lj(i)%sigma+lj(j)%sigma)
  157. epsmix=dsqrt(lj(i)%epslon*lj(j)%epslon)
  158.  
  159.  
  160. return
  161. end subroutine mixpar
  162.  
  163. subroutine num_der(valor,nu)
  164. ! Calcula las derivadas numéricas
  165.  
  166. !real(rk) :: valor(3,n_at), deltax, junk, r
  167. real(rk) :: valor, deltax, junk, r
  168. real(rk), dimension(2) :: pot
  169. integer, intent(in) :: nu
  170. integer :: i, j, k
  171.  
  172. !deltax = 1E-4_rk
  173. deltax = 0.0001_rk
  174.  
  175. write(*,*) coor
  176.  
  177. r = 0._rk
  178. do i = 1, 1
  179. do j = 1, 3
  180. coor(j, i) = coor(j, i) + deltax
  181. end do! j = 1, 3
  182.  
  183. do j = 1, 3
  184. write(*,*) "r(step) = ", r
  185. write(*,*) "coortest ", (coor(j,1) - coor(j,2))**2
  186. r = r + (coor(j,1) - coor(j,2))**2
  187. end do! k =
  188.  
  189. r = dsqrt(r)
  190. end do! i = 1, n_at
  191. write(*,*) "r = ", r
  192. pot(1) = 1._rk/r**nu
  193.  
  194. r = 0._rk
  195. do i = 1,n_at
  196. do j = 1, 3
  197. coor(j, i) = coor(j, i) -2._rk*deltax
  198. end do! j = 1, 3
  199.  
  200. do j = 1, 3
  201. r = r + (coor(j,1) - coor(j,2))**2
  202. end do! k =
  203.  
  204. r = dsqrt(r)
  205. end do! i = 1, n_at
  206. write(*,*) "r = ", r
  207. pot(2) = 1._rk/r**nu
  208.  
  209. !valor = (pot(1) - pot(2))/ (2._rk*deltax)
  210.  
  211. end subroutine !num_der
  212.  
  213. subroutine numderv2
  214. real(rk) :: forcen(3,n_at), p_ant(3,n_at), p_desp(3,n_at), dif(3,n_at), dist(3), ndev(3,n_at)
  215. real(rk), parameter :: deltax=1E-4_rk
  216. integer :: i , j, k, nu
  217. nu=2
  218. ndev=0._rk
  219. p_ant=0._rk
  220. p_desp=0._rk
  221. dif=0._rk
  222. do i=1, n_at
  223. !do i=1,3
  224. !x(j:i)=x(j:i)+deltax
  225. do k=i+1, n_at-1
  226. ! dist=rsq(xc(:,i),xc(:,j),boxlen)
  227. ! call rsq_cal(xc(:,i),xc(:,j),boxlen,dist,rvec,ndev)
  228. dist=xc(:,i)-xc(:,k)
  229. do j=1,3
  230. dist(j)=dist(j)-boxlen(j)*anint(dist(j)/boxlen(j))
  231. enddo!j
  232. dist=dot_product(dist,dist)
  233. dist=dsqrt(dist)
  234. do j=1,3
  235. p_ant(j,i)=p_ant(j,i)+1._rk/dist(j)**nu
  236. enddo!j
  237. enddo!k
  238. enddo!i
  239.  
  240. write(*,*) "derivadas", p_ant(:,1)
  241. stop
  242. endsubroutine !numder
  243.  
  244. end module force_cal1
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement