diff --git a/fields.f90 b/fields.f90 index 60f4a38..d0c84de 100644 --- a/fields.f90 +++ b/fields.f90 @@ -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 @@ -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' ! @@ -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") ! @@ -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,:) @@ -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 ! diff --git a/mol_xy2.f90 b/mol_xy2.f90 index 132afdb..cc1ce0f 100644 --- a/mol_xy2.f90 +++ b/mol_xy2.f90 @@ -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) @@ -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') @@ -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 @@ -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') @@ -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 @@ -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 @@ -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 @@ -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') ! @@ -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 @@ -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 ! @@ -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 diff --git a/mol_xy3.f90 b/mol_xy3.f90 index dc0821b..f2f7b02 100644 --- a/mol_xy3.f90 +++ b/mol_xy3.f90 @@ -35,9 +35,9 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) integer(ik) :: i,n,ix,jx,ioper character(cl) :: method ! - real(ark) :: rho_ark + real(ark) :: rho_ark,Inert0(3),resid(3,3) ! - integer(ik) :: ijk(3,2) + integer(ik) :: ijk(3,2),i0,j0,ipar,istep ! if (verbose>=4) write(out,"('ML_b0_XY3/start')") @@ -99,21 +99,23 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) ! end select ! - cosr = cos(rho) - sinr = sin(rho) + call generate_structure(rho,b0(1:4,1:3,0)) ! - b0(2,1,0) = re14*sinr - b0(2,2,0) = 0 - b0(2,3,0) = mX*re14*cosr/(Mtotal+mX) - b0(3,1,0) = -re14*sinr/2.0_ark - b0(3,2,0) = sqrt(3.0_ark)*re14*sinr/2.0_ark - b0(3,3,0) = mX*re14*cosr/(Mtotal+mX) - b0(4,1,0) = -re14*sinr/2.0_ark - b0(4,2,0) = -sqrt(3.0_ark)*re14*sinr/2.0_ark - b0(4,3,0) = mX*re14*cosr/(Mtotal+mX) - b0(1,1,0) = 0 - b0(1,2,0) = 0 - b0(1,3,0) = -Mtotal*re14*cosr/(Mtotal+mX) + !cosr = cos(rho) + !sinr = sin(rho) + ! + !b0(2,1,0) = re14*sinr + !b0(2,2,0) = 0 + !b0(2,3,0) = mX*re14*cosr/(Mtotal+mX) + !b0(3,1,0) = -re14*sinr/2.0_ark + !b0(3,2,0) = sqrt(3.0_ark)*re14*sinr/2.0_ark + !b0(3,3,0) = mX*re14*cosr/(Mtotal+mX) + !b0(4,1,0) = -re14*sinr/2.0_ark + !b0(4,2,0) = -sqrt(3.0_ark)*re14*sinr/2.0_ark + !b0(4,3,0) = mX*re14*cosr/(Mtotal+mX) + !b0(1,1,0) = 0 + !b0(1,2,0) = 0 + !b0(1,3,0) = -Mtotal*re14*cosr/(Mtotal+mX) ! ! Find center of mass ! @@ -149,16 +151,6 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) ! endif ! - if (verbose>=6) then - write(out,"(i6)") molec%natoms - ! - write(out,"(/'Bi',3x,3f14.8)") b0(1,:,i) - write(out,"( 'H',3x,3f14.8)") b0(2,:,i) - write(out,"( 'D',3x,3f14.8)") b0(3,:,i) - write(out,"( 'D',3x,3f14.8)") b0(4,:,i) - ! - endif - ! ! We can also simulate the reference structure at each point of rho, if npoints presented ! if (Npoints/=0) then @@ -180,90 +172,140 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) ! ! define the rho-type coordinate ! - ! transform0 = 0 ! forall (ix=1:3) transform0(ix,ix) = 1.0_ark ! - do i = 0,npoints - ! - rho = rho0+rho_i(i) - ! - re14 = molec%req(1) - ! - select case(trim(molec%coords_transform)) - ! - case('R-S-DELTA-MEP','R-PHI-DELTA-MEP','R-SYMPHI-DELTA-MEP') - ! - select case(trim(molec%meptype)) - case default - write (out,"('ML_b0_XY3: MEP type ',a,' unknown')") trim(molec%meptype) - stop 'ML_b0_XY3 - bad MEP type' - ! - case('MEP_H3OP') - ! - re14 = ML_mep_oh3p(rho) - ! - end select - ! - end select - ! - if (trim(molec%coords_transform)=='R-A2-A3-TAU') then - ! - phi = rho_i(i) - ! - hx = -sqrt(1.0_ark-2.0_ark*cos(phi) )/(-1.0_ark+cos(phi)) - hy = -cos(phi)/(-1.0_ark+cos(phi)) - ! - alpha = atan2(hx,hy) - ! - sinrho = 2.0_ark*sin(alpha*0.5_ark)/sqrt(3.0_ark) - ! - if ( sinrho>=1.0_ark ) then - rho = 0.5_ark*pi - else - rho = asin(sinrho) - endif + call generate_structure(rho0,a0(1:4,1:3)) + ! + ! symmetric XY3 case + ! + if (all(molec%AtomMasses(2:4)==mH1)) then + ! + do i = 0,npoints ! - if (phi>pi) rho = pi - rho + rho = rho0+rho_i(i) ! - endif - ! - cosr = cos(rho) - sinr = sin(rho) - ! - b0(2,1,i) = re14*sinr - b0(2,2,i) = 0.0_ark - b0(2,3,i) = mX*re14*cosr/(Mtotal+mX) - b0(3,1,i) = -re14*sinr/2.0_ark - b0(3,2,i) = sqrt(3.0_ark)*re14*sinr/2.0_ark - b0(3,3,i) = mX*re14*cosr/(Mtotal+mX) - b0(4,1,i) = -re14*sinr/2.0_ark - b0(4,2,i) = -sqrt(3.0_ark)*re14*sinr/2.0_ark - b0(4,3,i) = mX*re14*cosr/(Mtotal+mX) - b0(1,1,i) = 0.0_ark - b0(1,2,i) = 0.0_ark - b0(1,3,i) = -Mtotal*re14*cosr/(Mtotal+mX) - ! - do n = 1,3 - CM_shift = sum(b0(:,n,i)*molec%AtomMasses(:))/sum(molec%AtomMasses(:)) - b0(:,n,i) = b0(:,n,i) - CM_shift - enddo - ! - !forall(ix=1:4) b0(ix,:,i) = matmul(transpose(transform),b0(ix,:,i)) - ! - if (any(molec%AtomMasses(2:4)/=mH1)) then + re14 = molec%req(1) ! - ijk(1,1) = 2 ; ijk(1,2) = 3 - ijk(2,1) = 3 ; ijk(2,2) = 1 - ijk(3,1) = 1 ; ijk(3,2) = 2 + select case(trim(molec%coords_transform)) ! - method = 'ULEN' + case('R-S-DELTA-MEP','R-PHI-DELTA-MEP','R-SYMPHI-DELTA-MEP') + ! + select case(trim(molec%meptype)) + case default + write (out,"('ML_b0_XY3: MEP type ',a,' unknown')") trim(molec%meptype) + stop 'ML_b0_XY3 - bad MEP type' + ! + case('MEP_H3OP') + ! + re14 = ML_mep_oh3p(rho) + ! + end select + ! + end select ! - a0 = b0(:,:,i) + call generate_structure(rho,b0(1:4,1:3,i)) ! - call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i),transform,method=method) - !loop_xyz : do + enddo + ! + else + ! + ! non symmetric XY3 case (isotopologues) + ! + i0 = npoints/2 + ! + re14 = molec%req(1) + ! + do ipar = 0,1 + ! + istep = (-1)**(ipar+2) + ! + i = npoints/2 + ! + rho = rho0+rho_i(i) + ! + call generate_structure(rho,b0(1:4,1:3,i)) + ! + a0(:,:) = b0(:,:,i) + ! + call MLorienting_a0(molec%Natoms,molec%AtomMasses,a0,transform0) + ! + Inert0(1) = sum(molec%AtomMasses(:)*( a0(:,2)**2+ a0(:,3)**2) ) + Inert0(2) = sum(molec%AtomMasses(:)*( a0(:,1)**2+ a0(:,3)**2) ) + Inert0(3) = sum(molec%AtomMasses(:)*( a0(:,1)**2+ a0(:,2)**2) ) + ! + !do ix = 1,molec%Natoms + ! a0(ix,:) = matmul(transpose(transform),b0(ix,:,i)) + !enddo + ! + !b0(:,:,i) = a0(:,:) + ! + do j0 = 1,npoints/2 + ! + i = i + istep + ! + rho = rho0+rho_i(i) + ! + call generate_structure(rho,b0(1:4,1:3,i)) + ! + do ix = 1,molec%Natoms + a0(ix,:) = matmul(transpose(transform0),b0(ix,:,i)) + enddo + ! + call MLorienting_a0(molec%Natoms,molec%AtomMasses,a0,transform1) + ! + transform = matmul(transform1,transform0) + ! + do ix = 1,molec%Natoms + a0(ix,:) = matmul(transpose(transform),b0(ix,:,i)) + enddo + ! + b0(:,:,i) = a0(:,:) + ! + do ix = 1,3 + do jx = 1,3 + resid(ix,jx) = sum(molec%AtomMasses(:)*b0(:,ix,i)*b0(:,jx,i) ) + enddo + ! + enddo + ! + Inert(1) = sum(molec%AtomMasses(:)*( b0(:,2,i)**2+ b0(:,3,i)**2) ) + Inert(2) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,3,i)**2) ) + Inert(3) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,2,i)**2) ) + ! + select case (trim(molec%coords_transform)) + ! + case ('R-A2-A3-TAU') + ! + phi = rho_i(i) + ! + hx = -sqrt(1.0_ark-2.0_ark*cos(phi) )/(-1.0_ark+cos(phi)) + hy = -cos(phi)/(-1.0_ark+cos(phi)) + ! + alpha = atan2(hx,hy) + ! + sinrho = 2.0_ark*sin(alpha*0.5_ark)/sqrt(3.0_ark) + ! + if ( sinrho>=1.0_ark ) then + rho = 0.5_ark*pi + else + rho = asin(sinrho) + endif + ! + if (phi>pi) rho = pi - rho + ! + end select + ! + ijk(1,1) = 2 ; ijk(1,2) = 3 + ijk(2,1) = 3 ; ijk(2,2) = 1 + ijk(3,1) = 1 ; ijk(3,2) = 2 + ! + method = 'ULEN' + ! + a0 = b0(:,:,i) + ! + !call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i),transform1,method=method) ! loop_xyz : do ix =1,3 ! @@ -291,73 +333,173 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) enddo enddo loop_xyz ! - !enddo loop_xyz - ! - !do ioper = 1,0 - ! ! - ! tmat = ML_sym_rotat(mod(ioper,6),transform) - ! ! - ! if (ioper>6) tmat = transpose(tmat) - ! ! - ! forall(ix=1:3) unit(ix) = dot_product(transform0(ix,:),tmat(ix,:)) - ! ! - ! unit = unit - 1.0_ark - ! ! - ! if ( all( abs( unit(:) )<1.0e-1 ) ) exit - ! ! - !enddo - ! - transform0 = matmul(tmat,transform) - ! - forall(ix=1:4) b0(ix,:,i) = matmul(transpose(tmat),b0(ix,:,i)) - ! - !Inert(1) = sum(molec%AtomMasses(:)*( b0(:,2,i)**2+ b0(:,3,i)**2) ) - !Inert(2) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,3,i)**2) ) - !Inert(3) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,2,i)**2) ) - ! - endif - ! - ! Second Eckart equation - ! - do ix = 1,3 - do jx = 1,3 - if (ix/=jx) then - ! - a_t = sum(molec%AtomMasses(:)*b0(:,ix,i)*b0(:,jx,i) ) - ! - if (abs(a_t)>100.0_ark*small_) then - write(out,"('ML_b0_XY3: b0 is not a solution of Eckart 2 for ix,jx =',2i4,d18.8)") ix,jx,a_t - stop 'ML_b0_XY3: b0 is not solution of Eckart2' - endif + transform0 = matmul(tmat,transform) + ! + forall(ix=1:4) b0(ix,:,i) = matmul(transpose(tmat),b0(ix,:,i)) + ! + enddo + ! + !Inert(1) = sum(molec%AtomMasses(:)*( b0(:,2,i)**2+ b0(:,3,i)**2) ) + !Inert(2) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,3,i)**2) ) + !Inert(3) = sum(molec%AtomMasses(:)*( b0(:,1,i)**2+ b0(:,2,i)**2) ) + ! + !write(out,"(f15.8,1x,3f15.8)") rho,Inert(1:3) + ! + enddo + ! + endif + ! + endif + ! + do i = 0,npoints + ! + ! Second Eckart equation + ! + do ix = 1,3 + do jx = 1,3 + if (ix/=jx) then + ! + a_t = sum(molec%AtomMasses(:)*b0(:,ix,i)*b0(:,jx,i) ) + ! + if (abs(a_t)>100.0_ark*small_) then + write(out,"('ML_b0_XY3: b0 is not a solution of Eckart 2 for ix,jx =',2i4,d18.8)") ix,jx,a_t + stop 'ML_b0_XY3: b0 is not solution of Eckart2' endif - enddo - ! + endif enddo ! - if (verbose>=3) then - write(out,"(i6)") molec%natoms - ! - write(out,"(/'N',3x,3f14.8)") b0(1,:,i) - write(out,"( 'D',3x,3f14.8)") b0(2,:,i) - write(out,"( 'H',3x,3f14.8)") b0(3,:,i) - write(out,"( 'H',3x,3f14.8)") b0(4,:,i) - ! - endif - ! - !call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i),transform) - ! - !do n = 1,molec%Natoms - ! b0(n,:,i) = matmul(transpose(transform),b0(n,:,i)) - !enddo - ! - if (verbose>=5) then - write(out,"('b0',i4,12f12.8)") i,b0(:,:,i) - endif - ! enddo - endif + ! + if (verbose>=3) then + ! + write(out,"(i6)") molec%natoms + ! + 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) + write(out,"( a,3x,3f14.8)") trim(molec%zmatrix(4)%name),b0(4,:,i) + ! + endif + ! + if (verbose>=5) then + write(out,"('b0',i4,12f12.8)") i,b0(:,:,i) + endif + ! + enddo ! if (verbose>=4) write(out,"('ML_b0_XY3/end')") + ! + contains + ! + subroutine generate_structure(rho,b0) + ! + real(ark),intent(in) :: rho + real(ark),intent(out) :: b0(4,3) + real(ark) :: sinr,cosr,CM_shift,theta + ! + select case(trim(molec%coords_transform)) + ! + case('R-THETA-TAU') + ! + theta = molec%alphaeq(2) + ! + b0(1,1) = 0.0_ark + b0(1,2) = 0.0_ark + b0(1,3) = 0.0_ark + ! + b0(2,1) = 0.0_ark + b0(2,2) = 0.0_ark + b0(2,3) = molec%req(1) + ! + b0(3,1) = molec%req(2)*sin(theta)*cos(rho*0.5_ark) + b0(3,2) =-molec%req(2)*sin(theta)*sin(rho*0.5_ark) + b0(3,3) = molec%req(2)*cos(theta) + ! + b0(4,1) = molec%req(2)*sin(theta)*cos(rho*0.5_ark) + b0(4,2) = molec%req(2)*sin(theta)*sin(rho*0.5_ark) + b0(4,3) = molec%req(2)*cos(theta) + ! + case default + ! + cosr = cos(rho) + sinr = sin(rho) + ! + if (all(molec%AtomMasses(2:4)==mH1)) then + ! + b0(2,1) = re14*sinr + b0(2,2) = 0.0_ark + b0(2,3) = mX*re14*cosr/(Mtotal+mX) + b0(3,1) = -re14*sinr/2.0_ark + b0(3,2) = sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(3,3) = mX*re14*cosr/(Mtotal+mX) + b0(4,1) = -re14*sinr/2.0_ark + b0(4,2) = -sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(4,3) = mX*re14*cosr/(Mtotal+mX) + b0(1,1) = 0.0_ark + b0(1,2) = 0.0_ark + b0(1,3) = -Mtotal*re14*cosr/(Mtotal+mX) + ! + else + ! + select case(trim(molec%frame)) + ! + case default ! R1-Z-R2-X-Y + ! + b0(2,3) = re14*sinr + b0(2,1) = 0.0_ark + b0(2,2) = mX*re14*cosr/(Mtotal+mX) + b0(3,3) = -re14*sinr/2.0_ark + b0(3,1) = sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(3,2) = mX*re14*cosr/(Mtotal+mX) + b0(4,3) = -re14*sinr/2.0_ark + b0(4,1) = -sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(4,2) = mX*re14*cosr/(Mtotal+mX) + b0(1,3) = 0.0_ark + b0(1,1) = 0.0_ark + b0(1,2) = -Mtotal*re14*cosr/(Mtotal+mX) + ! + case ('R1-X-R2-Y-Z') + ! + b0(2,1) = re14*sinr + b0(2,2) = 0.0_ark + b0(2,3) = mX*re14*cosr/(Mtotal+mX) + b0(3,1) = -re14*sinr/2.0_ark + b0(3,2) = sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(3,3) = mX*re14*cosr/(Mtotal+mX) + b0(4,1) = -re14*sinr/2.0_ark + b0(4,2) = -sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(4,3) = mX*re14*cosr/(Mtotal+mX) + b0(1,1) = 0.0_ark + b0(1,2) = 0.0_ark + b0(1,3) = -Mtotal*re14*cosr/(Mtotal+mX) + ! + case ('R1-X-R2-Z-Y') + ! + b0(2,1) = re14*sinr + b0(2,3) = 0.0_ark + b0(2,2) =-mX*re14*cosr/(Mtotal+mX) + b0(3,1) = -re14*sinr/2.0_ark + b0(3,3) = sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(3,2) =-mX*re14*cosr/(Mtotal+mX) + b0(4,1) = -re14*sinr/2.0_ark + b0(4,3) = -sqrt(3.0_ark)*re14*sinr/2.0_ark + b0(4,2) =-mX*re14*cosr/(Mtotal+mX) + b0(1,1) = 0.0_ark + b0(1,3) = 0.0_ark + b0(1,2) = Mtotal*re14*cosr/(Mtotal+mX) + ! + end select + ! + endif + ! + end select + ! + do n = 1,3 + CM_shift = sum(b0(:,n)*molec%AtomMasses(:))/sum(molec%AtomMasses(:)) + b0(:,n) = b0(:,n) - CM_shift + enddo + ! + end subroutine end subroutine ML_b0_XY3 @@ -3105,12 +3247,46 @@ subroutine ML_rotsymmetry_XY3(J,K,tau,gamma,ideg) ! case('C2V','C2V(M)') ! - gamma = 0 - ideg = 1 - if (mod(K+2,2)==0.and.tau==0) gamma = 1 !; return - if (mod(K+2,2)==0.and.tau==1) gamma = 2 !; return - if (mod(K+2,2)/=0.and.tau==0) gamma = 3 !; return - if (mod(K+2,2)/=0.and.tau==1) gamma = 4 !; return + ! + select case(trim(molec%frame)) + ! + case default + ! + gamma = 0 + ideg = 1 + ! + ! R1-X-R2-Z-Y + ! + ! for the planer configuration, x is along the symmetry axis, z is in the plane and y is orhogonal to the plane + ! + if (mod(K+2,2)==0.and.tau==0) gamma = 1 !; return + if (mod(K+2,2)==0.and.tau==1) gamma = 2 !; return + if (mod(K+2,2)/=0.and.tau==0) gamma = 3 !; return + if (mod(K+2,2)/=0.and.tau==1) gamma = 4 !; return + ! + case('R1-X-R2-Y-Z') + ! + ! for the planer configuration, x is along the symmetry axis, y is in the plane and y is orhogonal to the plane + ! + gamma = 0 + ideg = 1 + if (mod(K+2,2)==0.and.tau==0) gamma = 1 !; return + if (mod(K+2,2)==0.and.tau==1) gamma = 2 !; return + if (mod(K+2,2)/=0.and.tau==0) gamma = 3 !; return + if (mod(K+2,2)/=0.and.tau==1) gamma = 4 !; return + ! + case('R1-Z-R2-X-Y') + ! + ! for the planer configuration, z is along the symmetry axis, x is in the plane and y is orhogonal to the plane + ! + gamma = 0 + ideg = 1 + if (mod(K+2,2)==0.and.tau==0) gamma = 1 !; return + if (mod(K+2,2)==0.and.tau==1) gamma = 2 !; return + if (mod(K+2,2)/=0.and.tau==0) gamma = 4 !; return + if (mod(K+2,2)/=0.and.tau==1) gamma = 3 !; return + ! + end select ! case('C','C(M)') ! diff --git a/mol_zxy2.f90 b/mol_zxy2.f90 index 2bba878..7f9c505 100644 --- a/mol_zxy2.f90 +++ b/mol_zxy2.f90 @@ -447,11 +447,14 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) real(ark), intent(in),optional :: rho_i(0:Npoints) real(ark), intent(out),optional :: rho_ref real(ark), intent(in),optional :: rho_borders(2) ! rhomim, rhomax - borders - - real(ark) :: a0(molec%Natoms,3),CM_shift,tau,cosalpha_2,costau,delta,cosrho,req(1:3),alphaeq(1:2),tau_ - real(ark) :: xieq(6),rho,theta,sint_2,theta12,beta,m1,m2,m3,m4,ax,ay - integer(ik) :: Nbonds,i,n - + ! + real(ark) :: a0(molec%Natoms,3),CM_shift,tau,cosalpha_2,costau,delta,cosrho,req(1:3),alphaeq(1:2),tau_ + real(ark) :: xieq(6),rho,theta,sint_2,theta12,beta,m1,m2,m3,m4,ax,ay + integer(ik) :: Nbonds,i,n + ! + integer(ik) :: ijk(3,2),ix,ioper + real(ark) :: phi,tmat(3,3),transform(3,3),transform1(3,3),transform0(3,3),unit(3) + character(cl) :: method if (verbose>=4) write(out,"('ML_b0_ZXY2/start')") @@ -535,6 +538,10 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) ! rho_ref = 0 ! pi ! + transform0 = 0 + ! + forall (ix=1:3) transform0(ix,ix) = 1.0_ark + ! do i = 0,npoints ! theta = molec%alphaeq(1) @@ -607,40 +614,88 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) req(1:3) = xieq(1:3) alphaeq(1:2) = xieq(1:2) ! - endif + endif ! if ( abs(req(1)-req(2))pi*5.0/3.0-small_) then + ! write(out,"('ML_b0_ZXY2 error: R-THETA-TAU, illegal tau+pi is outside (60..300)',f15.8)") tau*180.0/pi + ! stop 'ML_b0_ZXY2 error: R-THETA-TAU, illegal tau+pi is outside (60..300)' + !endif + ! + call generate_structure(theta,tau,b0(:,:,i)) + ! + !theta = acos( cos(tau)/(1.0_ark-cos(tau)) ) + ! + ijk(1,1) = 2 ; ijk(1,2) = 3 + ijk(2,1) = 3 ; ijk(2,2) = 1 + ijk(3,1) = 1 ; ijk(3,2) = 2 + ! + method = 'ULEN' + ! + a0 = b0(:,:,i) + ! + call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i),transform,method=method) + ! + loop_xyz : do ix =1,3 + ! + do ioper = 1,4 + ! + phi = real(ioper-1,ark)*pi*0.5_ark + ! + tmat(:,:) = 0 + ! + tmat(ix,ix) = 1.0_ark + ! + tmat(ijk(ix,1),ijk(ix,1)) = cos(phi) + tmat(ijk(ix,1),ijk(ix,2)) = sin(phi) + tmat(ijk(ix,2),ijk(ix,1)) =-sin(phi) + tmat(ijk(ix,2),ijk(ix,2)) = cos(phi) + ! + transform1 = matmul(tmat,transform) + ! + forall(ix=1:3) unit(ix) = dot_product(transform0(ix,:),transform1(ix,:)) + ! + unit = unit - 1.0_ark + ! + if ( all( abs( unit(:) )<1.0e-1 ) ) exit loop_xyz + ! + enddo + ! + enddo loop_xyz + ! + transform0 = matmul(tmat,transform) + ! + forall(ix=1:4) b0(ix,:,i) = matmul(transpose(tmat),b0(ix,:,i)) + ! + else + ! + b0(1,1,i) = 0.0_ark + b0(1,2,i) = 0.0_ark + b0(1,3,i) = 0.0_ark + ! + b0(2,1,i) = 0.0_ark + b0(2,2,i) = 0.0_ark + b0(2,3,i) = molec%req(1) + ! + b0(3,1,i) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) + b0(3,2,i) =-molec%req(2)*sin(theta)*sin(tau*0.5_ark) + b0(3,3,i) = molec%req(2)*cos(theta) + ! + b0(4,1,i) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) + b0(4,2,i) = molec%req(2)*sin(theta)*sin(tau*0.5_ark) + b0(4,3,i) = molec%req(2)*cos(theta) + ! + call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i)) + ! + ! Find center of mass + ! + do n = 1,3 + CM_shift = sum(b0(:,n,i)*molec%AtomMasses(:))/sum(molec%AtomMasses(:)) + b0(:,n,i) = b0(:,n,i) - CM_shift + enddo ! endif - - - ! - b0(1,1,i) = 0.0_ark - b0(1,2,i) = 0.0_ark - b0(1,3,i) = 0.0_ark - ! - b0(2,1,i) = 0.0_ark - b0(2,2,i) = 0.0_ark - b0(2,3,i) = molec%req(1) - ! - b0(3,1,i) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) - b0(3,2,i) =-molec%req(2)*sin(theta)*sin(tau*0.5_ark) - b0(3,3,i) = molec%req(2)*cos(theta) - ! - b0(4,1,i) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) - b0(4,2,i) = molec%req(2)*sin(theta)*sin(tau*0.5_ark) - b0(4,3,i) = molec%req(2)*cos(theta) - ! - call MLorienting_a0(molec%Natoms,molec%AtomMasses,b0(:,:,i)) - ! - ! Find center of mass - ! - do n = 1,3 - CM_shift = sum(b0(:,n,i)*molec%AtomMasses(:))/sum(molec%AtomMasses(:)) - b0(:,n,i) = b0(:,n,i) - CM_shift - enddo ! !do n = 1,molec%Natoms ! b0(n,:,i) = matmul(transpose(transform),b0(n,:,i)) @@ -665,6 +720,38 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders) endif ! if (verbose>=4) write(out,"('ML_b0_ZXY2/end')") + ! + contains + ! + subroutine generate_structure(theta,tau,b0) + + real(ark),intent(in) :: theta,tau + real(ark),intent(out) :: b0(4,3) + real(ark) :: sinr,cosr,CM_shift + ! + b0(1,1) = 0.0_ark + b0(1,2) = 0.0_ark + b0(1,3) = 0.0_ark + ! + b0(2,1) = 0.0_ark + b0(2,2) = 0.0_ark + b0(2,3) = molec%req(1) + ! + b0(3,1) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) + b0(3,2) =-molec%req(2)*sin(theta)*sin(tau*0.5_ark) + b0(3,3) = molec%req(2)*cos(theta) + ! + b0(4,1) = molec%req(2)*sin(theta)*cos(tau*0.5_ark) + b0(4,2) = molec%req(2)*sin(theta)*sin(tau*0.5_ark) + b0(4,3) = molec%req(2)*cos(theta) + ! + do n = 1,3 + CM_shift = sum(b0(:,n)*molec%AtomMasses(:))/sum(molec%AtomMasses(:)) + b0(:,n) = b0(:,n) - CM_shift + enddo + ! + end subroutine + end subroutine ML_b0_ZXY2 diff --git a/moltype.f90 b/moltype.f90 index e526aa3..af22ec6 100644 --- a/moltype.f90 +++ b/moltype.f90 @@ -155,6 +155,7 @@ subroutine calc_func(x, y) type(MLZmatrixT),pointer :: zmatrix(:) ! ! character(len=cl) :: coords_transform ! type of the coordinate transformation + character(len=cl) :: frame ! body fixed frame integer(ik),pointer :: dihedtype(:) ! dihedral type character(len=cl),pointer :: coordinates(:,:) ! Identifying the coordinate system, e.g. 'Cartesian', 'Bond-Angle' integer(ik) :: Nmodes,Ndihedrals,Natoms,Nbonds,Nangles,Ncoords @@ -315,7 +316,7 @@ subroutine MLinitialize_molec(Moltype,Coordinates,coords_transform,& AtomMasses,local_eq, & force_,forcename_,ifit_,pot_ind_,specparam,potentype,kinetic_type,& IO_primitive,chk_numerov_fname,& - symmetry_,rho_border,zmatrix_) + symmetry_,rho_border,zmatrix_,frame) character(len=cl),intent(in) :: Moltype @@ -337,6 +338,7 @@ subroutine MLinitialize_molec(Moltype,Coordinates,coords_transform,& character(len=cl),intent(in) :: IO_primitive,chk_numerov_fname real(ark) :: rho_border(2) ! rhomim, rhomax - borders type(MLZmatrixT),intent(in) :: zmatrix_(:) ! + character(len=cl),intent(in) :: frame ! coordinate transformation type ! integer(ik) :: alloc,Ncoords ! @@ -344,6 +346,7 @@ subroutine MLinitialize_molec(Moltype,Coordinates,coords_transform,& ! molec%moltype = Moltype molec%coords_transform = coords_transform + molec%frame = frame ! molec%Natoms = size(AtomMasses) molec%Ndihedrals = Ndihedrals diff --git a/pot_xy2.f90 b/pot_xy2.f90 index 14b1baf..d72e8c8 100644 --- a/pot_xy2.f90 +++ b/pot_xy2.f90 @@ -3981,7 +3981,7 @@ function MLloc2pqr_xy2(r) result(f) ! select case(trim(molec%coords_transform)) ! - case('R-RHO-Z','R-RHO-Z-M2-M3') + case('R-RHO-Z','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT') ! a0(2, 1) = -r(1) * cos(alpha_2) a0(2, 3) = -r(1) * sin(alpha_2) @@ -4109,7 +4109,7 @@ function MLloc2pqr_xyz(r) result(f) a0(3,3) = -(sin(v)*R2*sin(alpha)*mX+sin(v)*R2*sin(alpha)*mY-cos(v)*R2*cos(alpha)*mX-cos(v)*R2*cos(alpha)*mY+& cos(v)*mY*R1)/(mX+mY+mZ) ! - case('R-RHO-Z','R-RHO-Z-M2-M3') + case('R-RHO-Z','R-RHO-Z-M2-M3','R-RHO-Z-M2-M3-BISECT') ! alpha_2 = r(3)*0.5_ark !