Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- Module force_cal1
- use types
- use global_variables
- use math_subroutines
- use neigbours1
- use potpar1
- implicit none
- save
- Contains
- subroutine force_calc(v,xc,deriv)
- real(rk)::v,xc(3,n_at)
- real(rk)::epsmix,sigmix,sig2,sig6,sig12,vij,r1,fij(3)
- real(rk)::vcmix,dvcmix,dvccmix
- integer(ik)::i,j,k,ndev(3)
- logical::deriv
- if (pottyp.eq.0) then
- V=0.0_rk
- if (simtyp.eq.1) force=0.0_rk
- endif
- if (pottyp.eq.1) then
- rvec=0.0_rk
- v=0.0_rk
- force=0.0_rk
- P=0.0_rk
- der_boxlen=0.0_rk
- do i=1,n_at-1
- do j=i+1,n_at
- !r2=rsq(xc(:,i),xc(:,j),boxlen)
- call rsq_cal(xc(:,i),xc(:,j),boxlen,r2,rvec,ndev)
- v=r2**(-6)
- fij=2*(xc(:,i)-xc(:,j))/(r2**7)
- force(:,i)=force(:,i)+fij
- force(:,j)=force(:,j)-fij
- enddo
- enddo
- end if
- if (pottyp.eq.2) then
- call numderv2
- rvec=0.0_rk
- v=0.0_rk
- force=0.0_rk
- P=0.0_rk
- der_boxlen=0.0_rk
- do i=1,n_at-1
- do j=i+1,n_at
- r2=rsq(xc(:,i),xc(:,j),boxlen)
- v=r2**(-1)
- fij=2*(xc(:,i)-xc(:,j))/(r2**2)
- force(:,i)=force(:,i)+fij
- force(:,j)=force(:,j)-fij
- !derivadas
- write(*,*) "derivadas nuemricas"
- ! call numder(valor,2)
- write(*,*) "derivadas analitica"
- write(*,*) force(:,1)
- write(*,*) force(:,2)
- stop
- enddo
- enddo
- !potencial igual a 1 y 3 con esfera rigida y esfera blanda
- end if
- !Escribir potenciales y fuerzas para esfera blanda y rÃgida
- if (pottyp.eq.3) then
- rvec=0.0_rk
- v=0.0_rk
- force=0.0_rk
- P=0.0_rk
- der_boxlen=0.0_rk
- do i=1,n_at-1
- do j=i+1,n_at
- call rsq_cal(xc(:,i),xc(:,j),boxlen,r2,rvec,ndev)
- if (r2.gt.rc2) cycle
- r1=dsqrt(r2)
- if (att(i)%atnum.ne.att(j)%atnum) then
- call mixpar(i,j,sigmix,epsmix)
- vcmix=max(lj(i)%vc,lj(j)%vc)
- dvcmix=max(lj(i)%dvc,lj(j)%dvc)
- dvccmix=max(lj(i)%dvcc,lj(j)%dvcc)
- else
- epsmix=lj(i)%epslon
- sigmix=lj(i)%sigma
- vcmix=lj(i)%vc
- dvcmix=lj(i)%dvc
- dvccmix=lj(i)%dvcc
- endif
- sig2=sigmix**2/r2
- sig6=sig2*sig2*sig2
- sig12=sig6*sig6
- v=v+4.0_rk*epsmix*(sig12-sig6-vcmix+dvccmix*(r1-rcut))
- if (deriv) then
- fij=24.0_rk*epsmix*(2*sig12-sig6-dvcmix*r1)*rvec/r2
- P=P+dot_product(fij,rvec)
- force(:,i)=force(:,i)+fij
- force(:,j)=force(:,j)-fij
- do k=1,3
- if (ndev(k).ne.0) then
- der_boxlen(k)=der_boxlen(k)+ndev(k)*fij(k)
- endif
- enddo
- endif
- enddo
- enddo
- endif
- return
- end subroutine force_calc
- !subroutine graph
- !integer(ik)::i
- !real(rk)::sig2,sig6,sig12,v,rvec(3),r2,fij(3),r
- !do i=0,1200
- ! r=0.5_rk+0.01_rk*i
- ! r2=r*r
- ! rvec=(/0.0_rk,0.0_rk,r/)
- ! v=0.0_rk
- ! sig2=sigma**2/r2
- ! sig6=sig2*sig2*sig2
- ! sig12=sig6*sig6
- ! if (r2.le.rc2) then
- ! v=v+sig12-sig6-vc+dvcc*(r-rcut)
- ! fij=(2.0_rk*sig12-sig6-dvc*r)*rvec/r2
- ! v=4.0_rk*epslon*v
- ! fij=24.0_rk*epslon*v
- ! endif
- ! write(1334,*) r,V,fij(3)
- !enddo
- !end subroutine graph
- subroutine mixpar(i,j,sigmix,epsmix)
- integer(ik),intent(in)::i,j
- real(rk)::sigmix,epsmix
- sigmix=0.5_rk*(lj(i)%sigma+lj(j)%sigma)
- epsmix=dsqrt(lj(i)%epslon*lj(j)%epslon)
- return
- end subroutine mixpar
- subroutine num_der(valor,nu)
- ! Calcula las derivadas numéricas
- !real(rk) :: valor(3,n_at), deltax, junk, r
- real(rk) :: valor, deltax, junk, r
- real(rk), dimension(2) :: pot
- integer, intent(in) :: nu
- integer :: i, j, k
- !deltax = 1E-4_rk
- deltax = 0.0001_rk
- write(*,*) coor
- r = 0._rk
- do i = 1, 1
- do j = 1, 3
- coor(j, i) = coor(j, i) + deltax
- end do! j = 1, 3
- do j = 1, 3
- write(*,*) "r(step) = ", r
- write(*,*) "coortest ", (coor(j,1) - coor(j,2))**2
- r = r + (coor(j,1) - coor(j,2))**2
- end do! k =
- r = dsqrt(r)
- end do! i = 1, n_at
- write(*,*) "r = ", r
- pot(1) = 1._rk/r**nu
- r = 0._rk
- do i = 1,n_at
- do j = 1, 3
- coor(j, i) = coor(j, i) -2._rk*deltax
- end do! j = 1, 3
- do j = 1, 3
- r = r + (coor(j,1) - coor(j,2))**2
- end do! k =
- r = dsqrt(r)
- end do! i = 1, n_at
- write(*,*) "r = ", r
- pot(2) = 1._rk/r**nu
- !valor = (pot(1) - pot(2))/ (2._rk*deltax)
- end subroutine !num_der
- subroutine numderv2
- 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)
- real(rk), parameter :: deltax=1E-4_rk
- integer :: i , j, k, nu
- nu=2
- ndev=0._rk
- p_ant=0._rk
- p_desp=0._rk
- dif=0._rk
- do i=1, n_at
- !do i=1,3
- !x(j:i)=x(j:i)+deltax
- do k=i+1, n_at-1
- ! dist=rsq(xc(:,i),xc(:,j),boxlen)
- ! call rsq_cal(xc(:,i),xc(:,j),boxlen,dist,rvec,ndev)
- dist=xc(:,i)-xc(:,k)
- do j=1,3
- dist(j)=dist(j)-boxlen(j)*anint(dist(j)/boxlen(j))
- enddo!j
- dist=dot_product(dist,dist)
- dist=dsqrt(dist)
- do j=1,3
- p_ant(j,i)=p_ant(j,i)+1._rk/dist(j)**nu
- enddo!j
- enddo!k
- enddo!i
- write(*,*) "derivadas", p_ant(:,1)
- stop
- endsubroutine !numder
- end module force_cal1
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement