From b08161f65da613cf9f1c785e13e1bb09f58877a2 Mon Sep 17 00:00:00 2001 From: Sergey Yurchenko Date: Sun, 26 May 2024 22:35:02 +0100 Subject: [PATCH] CS2 PES by Huang implemented as a general type --- fields.f90 | 2 +- pot_CS2_ames1.f90 | 166 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 167 insertions(+), 1 deletion(-) create mode 100644 pot_CS2_ames1.f90 diff --git a/fields.f90 b/fields.f90 index 927df40..1feb84a 100644 --- a/fields.f90 +++ b/fields.f90 @@ -3594,7 +3594,7 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot) ! pot_form_compact = .true. ! - case("NAME") + case("NAME","POT_NAME","POT-NAME") ! pot_form_compact = .true. ! diff --git a/pot_CS2_ames1.f90 b/pot_CS2_ames1.f90 new file mode 100644 index 0000000..ef59573 --- /dev/null +++ b/pot_CS2_ames1.f90 @@ -0,0 +1,166 @@ +! +! This unit is for a user defined potential of CH3OH from Bowman (2007) +! +module pot_user + use accuracy + use moltype + + implicit none + + public MLdipole,MLpoten,ML_MEP,MLpoten_name + + private + + integer(ik), parameter :: verbose = 4 ! Verbosity level + ! + contains + ! + ! + function ML_MEP(dim,rho) result(f) + + integer(ik),intent(in) :: dim + real(ark),intent(in) :: rho + real(ark) :: f(dim) + ! + if (dim/=3) stop 'Illegal size of the function - must be 3' + ! + f(:) = molec%local_eq(:) + f(molec%Ncoords) = rho + + end function ML_MEP + ! + ! + recursive subroutine MLdipole(rank,ncoords,natoms,local,xyz,f) + ! + integer(ik),intent(in) :: rank,ncoords,natoms + real(ark),intent(in) :: local(ncoords),xyz(natoms,3) + real(ark),intent(out) :: f(rank) + ! + !f = MLdipole_xy4_dF(xyz) + ! + end subroutine MLdipole + ! + ! Check the potential name + subroutine MLpoten_name(name) + ! + character(len=cl),intent(in) :: name + character(len=cl),parameter :: poten_name = 'CS2_AMES1' + ! + if (poten_name/=trim(name)) then + write(out,"('a,a,a,a')") 'Wrong Potential ',trim(name),'; should be ',trim(poten_name) + endif + ! + write(out,"('a')") ' Using USER-tpye PES ',trim(poten_name) + ! + ! + end subroutine MLpoten_name + ! + ! Defining potential energy function) + function MLpoten(ncoords,natoms,local,xyz,force) result(f) + ! + integer(ik),intent(in) :: ncoords,natoms + real(ark),intent(in) :: local(ncoords) + real(ark),intent(in) :: xyz(natoms,3) + real(ark),intent(in) :: force(:) + real(ark) :: f + ! + f = MLpoten_cs2_ames1(ncoords,natoms,local,xyz,force) + ! + end function MLpoten + + + ! + function MLpoten_cs2_ames1(ncoords,natoms,local,xyz,force) result(f) + ! + implicit none + ! + integer(ik),intent(in) :: ncoords,natoms + real(ark),intent(in) :: local(ncoords) + real(ark),intent(in) :: xyz(natoms,3) + real(ark),intent(in) :: force(:) + real(ark) :: f + ! + integer(ik) :: i,k,Ncoe + real(ark) :: alpha1,alpha2,de1,de2,rcs1,rcs2,rcsref,rrefa,rrefb,rcsd1,rcsd2,ang2 + ! + real(ark) :: emin,etmp2,sumx0,ax,enetmp2,etmp1,enetmp3 + ! + real(ark) :: str1,str2,enetmp1,atp1,edamp11,edamp12,edamp1,V + real(ark) :: v0,r1,r2,a3,xx1,xx2,ang,sumstr2,sumstr4,angref,angx,bdamp2,bdamp4 + real(ark) :: edp1,edp2,edp3,edp4,edp5,edp6,De_1,De_2,De3,De_3 + ! + real(ark) :: v1(3),v2(3),x1,x2,cosalpha1 + ! + rcsref = force(1) + alpha1 = force(2) + alpha2 = alpha1 + De1 = force(3) + De2 = De1 + De_1 = force(4) + De_2 = De_1 + De3 = force(5) + De_3 = force(6) + edp1 = force(7) + edp2 = force(8) + edp3 = force(9) + edp4 = force(10) + edp5 = force(11) + edp6 = force(12) + Emin = force(13)*219474.63067_ark + ! + rcs1 = local(1) + rcs2 = local(2) + ! + ang = local(3) + ! + r1=rcs1-rcsref + r2=rcs2-rcsref + ! + a3=1.0_ark+cos(ang) + ! + Ncoe = molec%parmax + ! + v0=0 + do i=14,Ncoe + ! + v0=v0+force(i)*r1**molec%pot_ind(1,i)*r2**molec%pot_ind(2,i)*a3**(molec%pot_ind(3,i)) + if(molec%pot_ind(1,i).ne.molec%pot_ind(2,i))then + v0=v0+force(i)*r2**molec%pot_ind(1,i)*r1**molec%pot_ind(2,i)*a3**(molec%pot_ind(3,i)) + end if + ! + end do + ! + rrefa=rcsref; rrefb=rcsref + rcsd1=rcs1-rrefa; rcsd2=rcs2-rrefb + rcsd1=rcsd1*rcsd1; rcsd2=rcsd2*rcsd2 + sumstr2=rcsd1+rcsd2 + sumstr4=rcsd1*rcsd1+rcsd2*rcsd2 + ! + r1=rcs1-rcsref + r2=rcs2-rcsref + ! + ang2=pi-ang + angref=150.d0*pi/180.d0-pi + angx=min(-abs(ang2)-angref,0.0_ark) + bdamp2=angx**2; bdamp4=angx**4 + etmp2=exp(edp1*sumstr2+edp2*sumstr4+edp3*bdamp2+edp4*bdamp4) + ! + V0=V0*etmp2 + ! + enetmp1=De1*(1.0_ark-exp(-alpha1*(rcs1-rcsref)))**2 + De_1*(1.0_ark-exp(-alpha1*(rcs1-rcsref)))**4 + enetmp2=De2*(1.0_ark-exp(-alpha2*(rcs2-rcsref)))**2 + De_2*(1.0_ark-exp(-alpha2*(rcs2-rcsref)))**4 + ! + ax=(pi-ang)*0.5_ark + enetmp3=De3*sin(ax)**2 + De_3*sin(ax)**4 !bending simulation + edamp1=exp(edp5*sumstr2+1.d0*edp6*sumstr4) + enetmp3=edamp1*enetmp3 + ! + etmp1=enetmp1+enetmp2+enetmp3 + ! + V=V0-emin+etmp1 + ! + f = V + ! + end function MLpoten_cs2_ames1 + +end module pot_user