Skip to content

Commit

Permalink
frame is introduced; the ND2H and NH2D cases should be now working wi…
Browse files Browse the repository at this point in the history 
…th the standard XY3 R-delta-theta case

atomic names are printed out automatically
  • Loading branch information
@Trovemaster
Trovemaster committed Sep 29, 2023
1 parent bc1cfa5 commit fb25db8
Show file tree
Hide file tree
Showing 6 changed files with 530 additions and 220 deletions.
18 changes: 14 additions & 4 deletions fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,8 @@ module fields
character(len=cl) :: Moltype ! Identifying type of the Molecule (e.g. XY3)
character(len=cl),pointer :: Coordinates(:,:) ! Identifying the coordinate system, e.g. 'normal'; kinetic and potential part can be different
character(len=cl) :: internal_coords ! type of the internal coordinates: linear,harmonic,local
character(len=cl) :: coords_transform ! type of the coordinate transformation: linear,morse-xi-s-rho,...
character(len=cl) :: coords_transform = 'LINEAR' ! type of the coordinate transformation: linear,morse-xi-s-rho,...
character(len=cl) :: frame = "DEFAULT" ! body fixed frame
character(len=cl) :: symmetry ! molecular symmetry
integer(ik):: Natoms ! Number of atoms
integer(ik):: Nmodes ! Number of modes = 3*Natoms-6
Expand Down Expand Up @@ -724,7 +725,6 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
job%vib_contract = .false.
!
trove%internal_coords = 'LINEARIZED'
trove%coords_transform = 'LINEAR'
job%PTtype = 'POWERS'
if (manifold==1) job%PTtype = 'DIAGONAL'
!
Expand Down Expand Up @@ -1659,7 +1659,15 @@ subroutine FLReadInput(NPTorder,Npolyads,Natoms,Nmodes,Jrot)
!
call readu(w)
!
trove%coords_transform = trim(w)
trove%coords_transform = trim(w)
!
if (trim(trove%frame)=='DEFAULT') trove%frame = trim(w)
!
case("FRAME")
!
call readu(w)
!
trove%frame = trim(w)
!
case("MOLTYPE")
!
Expand Down Expand Up @@ -5302,6 +5310,8 @@ subroutine FLsetMolecule
!
write(out,"('The coordinate-transformation type (coords_transform):',a)") trim(trove%coords_transform)
!
write(out,"('The body-fixed frame is ',a)") trim(trove%frame)
!
write(out,"('Bonds connections are ')")
do ibond = 1,trove%Nbonds
write(out,"(i4,' - ',2i4)") ibond,trove%bonds(ibond,:)
Expand Down Expand Up @@ -5497,7 +5507,7 @@ subroutine FLsetMolecule
trove%mass,trove%local_eq,&
force,forcename,ifit,pot_ind,trove%specparam,trove%potentype,trove%kinetic_type,&
trove%IO_primitive,trove%chk_numerov_fname,&
trove%symmetry,trove%rho_border,trove%zmatrix)
trove%symmetry,trove%rho_border,trove%zmatrix,trove%frame)
!
! define the potential function method
!
Expand Down
58 changes: 46 additions & 12 deletions mol_xy2.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ function ML_coordinate_transform_XY2(src,ndst,direct) result (dst)
dst(:) = dst(:) + molec%local_eq(:)
endif
!
case('R-RHO','R-RHO-Z','R-RHO-Z-ECKART','R-RHO-Z-M2-M3')
case('R-RHO','R-RHO-Z','R-RHO-Z-ECKART','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT')
!
if (direct) then
dst(1:2) = dsrc(1:2)
Expand Down Expand Up @@ -718,7 +718,7 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
!b0(2,1,0) = (m1+m3)/m*re23*cos(e)+m3/m*re13*cos(rho-e)


select case(trim(molec%coords_transform))
select case(trim(molec%frame))
case default
!
case('R-RHO-Z','R-PHI-RHO-Z','R-RHO-Z-ECKART','R-RHO-HALF')
Expand Down Expand Up @@ -764,6 +764,26 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
b0(3,2,0) = 0.0_ark
b0(3,1,0) =-m1/m*re13*cos(alphae_h)
!
case('R-RHO-Z-M2-M3-BISECT')
!
if (Nangles==0) then
stop 'R-RHO-Z-M2-M3-BISECT does not work for Nangles = 0'
endif
!
alphae_h = alphaeq*0.5_ark
!
b0(1,1,0) = cos(alphae_h)*(m2*re13+m3*re23)/(m1+m2+m3)
b0(1,2,0) = 0
b0(1,3,0) = sin(alphae_h)*(m2*re13-m3*re23)/(m1+m2+m3)
!
b0(2,1,0) =-cos(alphae_h)*(re13*m1+re13*m3-m3*re23)/(m1+m2+m3)
b0(2,2,0) = 0
b0(2,3,0) =-sin(alphae_h)*(re13*m1+re13*m3+m3*re23)/(m1+m2+m3)
!
b0(3,1,0) =-cos(alphae_h)*(re23*m1+re23*m2-m2*re13)/(m1+m2+m3)
b0(3,2,0) = 0
b0(3,3,0) = sin(alphae_h)*(re23*m1+re23*m2+m2*re13)/(m1+m2+m3)
!
case('R-R-R')
!
b0(1,2,0) = 0.0_ark
Expand Down Expand Up @@ -944,9 +964,9 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
! stop 'ML_b0_XY2: rho_borders or rho_ref not specified '
!endif
!
select case(trim(molec%coords_transform))
select case(trim(molec%frame))
case default
write (out,"('ML_b0_XY2: coord. type ',a,' unknown')") trim(molec%coords_transform)
write (out,"('ML_b0_XY2: coord. type ',a,' unknown')") trim(molec%frame)
stop 'ML_b0_XY2 - bad coord. type'
!
case('LINEAR')
Expand All @@ -962,7 +982,7 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
a02 = (m1/m)
!
case('R-RHO','R-EXPRHO','R-RHO-Z','R12-R','R12-RHO','R13-RHO','R-PHI-RHO','R-PHI-RHO-Z','R-PHI1-PHI2-Z','R-PHI1-Z',&
'R-S1-S2-Z','R-ALPHA-THETA-Z','R-RHO-Z-ECKART','R-RHO-Z-M2-M3')
'R-S1-S2-Z','R-ALPHA-THETA-Z','R-RHO-Z-ECKART','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT')
!
rho_ref = 0.0_ark
rho0 = 0.0_ark
Expand Down Expand Up @@ -1005,7 +1025,7 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
!
select case(trim(molec%coords_transform))
case ('R-RHO','R12-RHO','R13-RHO','R-RHO-Z','R-PHI-RHO','R-PHI-RHO-Z',&
'R1-Z-R2-RHO','R-RHO-Z-ECKART','R1-Z-R2-RHO-ECKART','RADAU-R-ALPHA-Z','R-RHO-Z-M2-M3')
'R1-Z-R2-RHO','R-RHO-Z-ECKART','R1-Z-R2-RHO-ECKART','RADAU-R-ALPHA-Z','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT')
alpha = pi-rho
case('R-RHO-HALF')
alpha = pi-rho*2.0_ark
Expand All @@ -1017,7 +1037,7 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
!
!if(trim(molec%coords_transform)=='R-EXPRHO') alpha = pi-(1.0_ark+rho**2)/(max(rho,small_))
!
if(trim(molec%coords_transform)=='RADAU') then
if(trim(molec%frame)=='RADAU') then
!
!
!f1= 1.0_ark/g1
Expand Down Expand Up @@ -1110,7 +1130,7 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
b0(3,2,i) = 0.0_ark
b0(3,3,i) =-m1/m*re13*cos(alphae_h)
!
select case(trim(molec%coords_transform))
select case(trim(molec%frame))
!
case('R-RHO-Z','R-PHI-RHO-Z','R-RHO-Z-ECKART','R-RHO-HALF')
!
Expand Down Expand Up @@ -1140,6 +1160,20 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
b0(3,2,i) = 0.0_ark
b0(3,3,i) = re13*sin(alphae_h)*(m1+m2+m2)/m
!
case('R-RHO-Z-M2-M3-BISECT')
!
b0(1,1,i) = cos(alphae_h)*(m2*re13+m3*re23)/(m1+m2+m3)
b0(1,2,i) = 0
b0(1,3,i) = sin(alphae_h)*(m2*re13-m3*re23)/(m1+m2+m3)
!
b0(2,1,i) =-cos(alphae_h)*(re13*m1+re13*m3-m3*re23)/(m1+m2+m3)
b0(2,2,i) = 0
b0(2,3,i) =-sin(alphae_h)*(re13*m1+re13*m3+m3*re23)/(m1+m2+m3)
!
b0(3,1,i) =-cos(alphae_h)*(re23*m1+re23*m2-m2*re13)/(m1+m2+m3)
b0(3,2,i) = 0
b0(3,3,i) = sin(alphae_h)*(re23*m1+re23*m2+m2*re13)/(m1+m2+m3)
!
case('R-R-R')
!
b0(1,2,i) = 0.0_ark
Expand Down Expand Up @@ -1332,9 +1366,9 @@ subroutine ML_b0_XY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
if (verbose>=4) then
write(out,"(i6)") molec%natoms
!
write(out,"(/'C',3x,3f14.8)") b0(1,:,i) !*bohr
write(out,"( 'O',3x,3f14.8)") b0(2,:,i) !*bohr
write(out,"( 'O',3x,3f14.8)") b0(3,:,i) !*bohr
write(out,"(/a,3x,3f14.8)") trim(molec%zmatrix(1)%name),b0(1,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(2)%name),b0(2,:,i)
write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(3)%name),b0(3,:,i)
!
endif
!
Expand Down Expand Up @@ -2758,7 +2792,7 @@ subroutine ML_symmetry_transformation_XY2(ioper,natoms,src,dst)
!
end select
!
case('R1-Z-R2-ALPHA','R1-R2-ALPHA-Z','R1-Z-R2-RHO','R1-Z-R2-RHO-ECKART','R-RHO-Z-M2-M3')
case('R1-Z-R2-ALPHA','R1-R2-ALPHA-Z','R1-Z-R2-RHO','R1-Z-R2-RHO-ECKART','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT')
!
select case(trim(molec%symmetry))
case default
Expand Down
Loading

0 comments on commit fb25db8

Please sign in to comment.