C Fortran-90 program to convert Cartesian cooordinates of two water C molecules to two sets of Euler angles and an intermolecular C distance. C C Program is for use with SAPT5s.f and assumes, just like SAPT5s.f, C that the standard water is in the xz-plane with the oxygen on the C positive z-axis and the center of mass in the origin. C C C Input: C - namelist/input/ (from standard I/O): C coord_file = name of file containing Cartesian C coordinates of water dimer (char*60). C Default: 'coordinates.dat' C C output_file = name of file containing output (is input file C to sapt5s.f). Default: 'example.dat' C C what_units = length unit (string: 'angstrom' or 'bohr') C Default: bohr C C Example namelist: (do not copy the exclamation marks) C------------------------------------------------------ ! &input coord_file = 'flop.dat', ! output_file = 'sample.dat', ! what_units = 'angstrom' / C------------------------------------------------------- C C - coord_file contains free format lines: C atom, x, y, z C ... C atom, x, y, z C C Example coord_file ('flop.dat'). Do not copy exclamation marks. C -------------------------------------------------------------- !* Coordinates obtained from F. van Duijneveldt: !* Monomer 1 ! O1 0.873683 0.000000 -1.246754 ! h1 0.287800 0.000000 -2.003704 ! H2 0.287800 0.000000 -0.489804 !* Monomer 2 ! o2 -0.778985 0.000000 1.131754 ! H3 -0.666591 0.756950 1.706755 ! h4 -0.666591 -0.756950 1.706754 C----------------------------------------------------------------- C C Notes on structure of coord_file (ASCII) C ---------------------------------------- C - Lines on coord_file beginning with '*' are skipped (allows for C comments). C - Empty lines are skipped C - Atom is a string starting with O, H, or D. C - First three non-comment lines --> first monomer C - Second three non-comment lines --> second monomer C - Specification of monomers must start with oxygen atoms. C C Output: C - Monomer geometries, Euler angles, distance on standard I/O. C - Distance, Euler angles on file with name 'output_file' C C Author: P.E.S. Wormer, February 2001. MODULE constants implicit none character*60, public :: coord_file character*60, public :: output_file integer*4, public, parameter :: maxatoms = 6 real*8, public :: units real*8, public, parameter :: pi = 3.1415926535897932d0 real*8, public, parameter :: angle_HOH =104.69d0 real*8, public, parameter :: OH_length = 1.836106337d0 ! bohr real*8, public, parameter :: toler_l = OH_length/10.d0 real*8, public, parameter :: toler_a = angle_HOH/10.d0 real*8, public, parameter :: bohr = 0.5291772083d0 real*8, public, parameter :: m_H = 1.007825d0 real*8, public, parameter :: m_D = 2.0140d0 real*8, public, parameter :: m_O = 15.99491d0 real*8, public, parameter :: smallest = tiny(units) real*8, public, parameter :: mac_prec = epsilon(units) end module constants MODULE vector_functions contains FUNCTION cart2polar(R) C From Cartesians to polars (radians). European definition of C polars. C Input 3-vector R. C Output vector (/ Length_R, theta, phi /). use constants, only: pi, smallest, mac_prec implicit none real*8 st, ct, sp, cp, . theta, phi, . Length_R, Length_Rp, Length_diff, . cart2polar(3), . R(3), . diff(3) intent(in) R Length_R = sqrt(dot_product(R,R)) if ( abs(Length_R) < smallest ) then theta = 0.d0 phi = 0.d0 return endif theta = acos( R(3)/Length_R ) Length_Rp = sqrt( dot_product( R(1:2),R(1:2) ) ) if (abs(Length_Rp) < smallest ) then phi = 0.d0 else phi = acos(R(1)/Length_Rp) if ( R(2) < -smallest ) then phi = 2.d0*pi - phi endif endif cart2polar = (/Length_R, theta, phi/) C--Check * st = sin(theta) * ct = cos(theta) * sp = sin(phi) * cp = cos(phi) * * diff = R - Length_R*(/st*cp, st*sp, ct/) * Length_diff = sqrt(dot_product(diff, diff)) * * if (Length_diff > mac_prec*1.d7 ) then * write(*,*) ' Error in cart2polar, length differences: ', * . Length_diff * stop 16 * endif END FUNCTION cart2polar FUNCTION euler_rotmat(alpha, beta, gamma) C Euler rotation matrix R(alpha, beta, gamma); active convention. implicit none real*8 alpha, beta, gamma, . euler_rotmat(3,3), . Z1(3,3), Y(3,3), Z(3,3) intent(in) alpha, beta, gamma Z = 0.d0; Y = 0.d0; Z1 = 0.d0; Z(1,1) = cos(gamma) Z(2,1) = sin(gamma) Z(1,2) = -sin(gamma) Z(2,2) = cos(gamma) Z(3,3) = 1.d0 Y(1,1) = cos(beta) Y(3,1) = -sin(beta) Y(2,2) = 1.d0 Y(1,3) = sin(beta) Y(3,3) = cos(beta) Z1(1,1) = cos(alpha) Z1(2,1) = sin(alpha) Z1(1,2) = -sin(alpha) Z1(2,2) = cos(alpha) Z1(3,3) = 1.d0 euler_rotmat = matmul(matmul(Z1,Y),Z) end function euler_rotmat FUNCTION cross(a,b) C Cross product a x b implicit none real*8 a(3), b(3), cross(3) intent(in) a, b cross(1) = a(2)*b(3) - a(3)*b(2) cross(2) = a(3)*b(1) - a(1)*b(3) cross(3) = a(1)*b(2) - a(2)*b(1) end function cross end module vector_functions PROGRAM CART2EULER use constants, only: maxatoms implicit none real*8 coords(3, maxatoms) character*4 atoms (maxatoms) integer natoms Call rdstdio ! namelist from standard i/o Call getcoords(coords, atoms, natoms)! coordinates from coord_file C--Overall rotation of system such that RAB coincides with lab- C fixed z-axis. Modify array coords. Print R_AB on standard I/O. Call rot_overall(coords, atoms, natoms) C--Get monomer geometries. Print Euler angles of monomers on standard C I/O. Call get_monomers(coords) END SUBROUTINE RDSTDIO C Read namelist from standard I/O use constants, only: coord_file, output_file, units, bohr implicit none character*8 what_units namelist/input/coord_file, output_file, what_units C--Namelist defaults coord_file = 'coordinates.dat' output_file = 'example.dat' what_units = 'bohr' read(*,input) write(*,'(t10, a/t10,a/ 3(2a/) )') . ' Namelist input ', . ' ============== ', . ' coord_file = ', trim(adjustl(coord_file)), . ' output_file = ', trim(adjustl(output_file)), . ' what_units = ', what_units C--Convert to bohr if (what_units(1:1) .eq. 'a' .or. what_units(1:1) .eq. 'A' ) . then units = 1.d0/bohr elseif (what_units(1:1) .eq. 'b' .or. what_units(1:1) .eq. 'B' ) . then units = 1.d0 else write(*,'(/a)') . ' Wrong length units, must be angstrom or bohr ' stop 16 endif END SUBROUTINE getcoords(coords, atoms, natoms) C Get Cartesian coordinates from 'coord_file'. C Return in array atom(char*4) and coords(3, natoms) use constants, only: maxatoms, coord_file, units implicit none integer natoms, i character*1, parameter :: comment = '*' real*8 coords(3, maxatoms) character*4 atoms(maxatoms) character*100 line intent(out) coords, atoms, natoms open(unit=1, file=coord_file,form='formatted') natoms = 0 10 read(unit= 1,fmt='(a100)',end=20) line if ( (line(1:1) .eq. comment) .or. (len(trim(line)) .eq. 0)) then go to 10 else if (natoms .eq. maxatoms) then write(*,'(//a/a)') ' *** Array COORDS too small. ', . ' Increase maxatoms and recompile.' stop 16 endif natoms = natoms + 1 read(line,*) atoms(natoms), coords(:,natoms) go to 10 endif 20 write(*,'(t15,a/t15,a)') . ' Input coordinates ', . ' ================= ' do i = 1,natoms write(*,'(1x, a, 3(f12.6,1x) )') . atoms(i), coords(:,i) enddo C--Check the atoms if (atoms(1)(1:1) .ne. 'o' .and. atoms(1)(1:1) .ne. 'O' ) then write(*,*) ' First atom must be oxygen ' stop 16 endif if (atoms(2)(1:1) .ne. 'h' .and. atoms(2)(1:1) .ne. 'H' . .and. atoms(2)(1:1) .ne. 'd' .and. atoms(2)(1:1) .ne. 'D' ) then write(*,*) ' Second atom must be proton or deuteron ' stop 16 endif if (atoms(3)(1:1) .ne. 'h' .and. atoms(3)(1:1) .ne. 'H' . .and. atoms(3)(1:1) .ne. 'd' .and. atoms(3)(1:1) .ne. 'D' ) then write(*,*) ' Third atom must be proton or deuteron ' stop 16 endif if (atoms(4)(1:1) .ne. 'o' .and. atoms(4)(1:1) .ne. 'O' ) then write(*,*) ' Fourth atom must be oxygen ' stop 16 endif if (atoms(5)(1:1) .ne. 'h' .and. atoms(5)(1:1) .ne. 'H' . .and. atoms(5)(1:1) .ne. 'd' .and. atoms(5)(1:1) .ne. 'D' ) then write(*,*) ' Fifth atom must be proton or deuteron ' stop 16 endif if (atoms(6)(1:1) .ne. 'h' .and. atoms(6)(1:1) .ne. 'H' . .and. atoms(6)(1:1) .ne. 'd' .and. atoms(6)(1:1) .ne. 'D' ) then write(*,*) ' Sixth atom must be proton or deuteron ' stop 16 endif C--To bohr: coords = coords*units END SUBROUTINE rot_overall(coords, atoms, natoms) C Overall rotation of dimer use constants, only: maxatoms, m_H, m_D, m_O, bohr, mac_prec, . output_file use vector_functions, only: cart2polar, euler_rotmat implicit none integer i, natoms real*8 m_tot_A, . m_tot_B, . coords (3, maxatoms), . masses (maxatoms), . com_A (3), . com_B (3), . RAB (3), . RABnew (3), . polarAB(3), . rotR (3,3) character*4 atoms(maxatoms) intent(in) atoms, natoms intent(inout) coords C--Set array masses do i = 1, natoms if ( (atoms(i)(1:1).eq.'h') .or. (atoms(i)(1:1).eq.'H')) . then masses(i) = m_H elseif ( (atoms(i)(1:1).eq.'d') .or. (atoms(i)(1:1).eq.'D')) . then masses(i) = m_D elseif ( (atoms(i)(1:1).eq.'o') .or. (atoms(i)(1:1).eq.'O')) . then masses(i) = m_O endif enddo C--Set centers of mass m_tot_A = sum(masses(1:3)) ! total mass of A m_tot_B = sum(masses(4:6)) ! total mass of B com_A = matmul(coords(:,1:3), masses(1:3))/m_tot_A com_B = matmul(coords(:,4:6), masses(4:6))/m_tot_B C--Matrix rotR rotates RAB such that it lies along the z-axis RAB = com_B - com_A ! points from A to B polarAB = cart2polar(RAB) ! returns R, theta, phi rotR = transpose(euler_rotmat( polarAB(3), polarAB(2), 0.d0)) RABnew = matmul(rotR, RAB) ! Rotate RAB coords = matmul(rotR, coords) ! Rotate total system C--Check if ( sqrt(dot_product(RABnew(1:2), RABnew(1:2))) > mac_prec*1.d3) . then write(*,*) ' R not along z-axis, RABnew = ', RABnew stop 16 endif write(*,'(/a,f14.8,a,2x,f14.8,a)') . ' Intermolecular distance: ', polarAB(1), ' [bohr]', . polarAB(1)*bohr, ' [Angstrom]' open(unit=2, file=output_file,form='formatted') write(2,'(i1)') 1 write(2,*) polarAB(1) END SUBROUTINE get_monomers(coords) C Get the geometry and Euler angles of the two monomers use constants, only: maxatoms, pi, OH_length, angle_HOH, . toler_a, toler_l use vector_functions, only: cross implicit none real*8 l_OH1, l_OH2, l_OH3, l_OH4, . angle1, angle2, . alphaA, betaA, gammaA, . alphaB, betaB, gammaB, . coords(3, maxatoms), . OH1 (3), . OH2 (3), . OH3 (3), . OH4 (3), . monA (3,3), . monB (3,3) integer i intent(in) coords C--Get the geometry of the monomers OH1 = (coords(:,2) - coords(:,1)) OH2 = (coords(:,3) - coords(:,1)) OH3 = (coords(:,5) - coords(:,4)) OH4 = (coords(:,6) - coords(:,4)) l_OH1 = sqrt(dot_product(OH1, OH1)) l_OH2 = sqrt(dot_product(OH2, OH2)) l_OH3 = sqrt(dot_product(OH3, OH3)) l_OH4 = sqrt(dot_product(OH4, OH4)) write(*,'(/a/a/a,a/ (a,f12.7,a,f12.7) )') . ' Monomer Geometries [bohr/degrees]', . ' =================================', . ' A ', ' B ', . ' OH1 = ',l_OH1, ' OH3 = ', l_OH3, . ' OH2 = ',l_OH2, ' OH4 = ', l_OH4 OH1 = OH1/l_OH1 OH2 = OH2/l_OH2 OH3 = OH3/l_OH3 OH4 = OH4/l_OH4 angle1 = acos(dot_product(OH1,OH2))*180.d0/pi angle2 = acos(dot_product(OH3,OH4))*180.d0/pi C--Continue writing ... write(*,'( a, f9.4, a, f9.4 )') . ' angle1 = ', angle1, ' angle2 = ', angle2 C--Error checks if ( abs(OH_length-l_OH1) > toler_l ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: O-H1 does not have the right length ', . ' must be: ', OH_length, ' is ', l_OH1 endif if ( abs(OH_length-l_OH2) > toler_l ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: O-H2 does not have the right length ', . ' must be: ', OH_length, ' is ', l_OH2 endif if ( abs(OH_length-l_OH3) > toler_l ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: O-H3 does not have the right length ', . ' must be: ', OH_length, ' is ', l_OH3 endif if ( abs(OH_length-l_OH4) > toler_l ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: O-H4 does not have the right length ', . ' must be: ', OH_length, ' is ', l_OH4 endif if ( abs(angle_HOH-angle1) > toler_a ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: First H-O-H angle not correct ', . ' must be: ',angle_HOH , ' is ', angle1 endif if ( abs(angle_HOH-angle2) > toler_a ) then write(*,'(/a/a,f12.7,a,f12.7)') . ' *** WARNING: Second H-O-H angle not correct ', . ' must be: ', angle_HOH , ' is ', angle2 endif monA(:,3) = -(OH1 + OH2) monA(:,1) = (OH1 - OH2) monA(:,3) = monA(:,3)/sqrt(dot_product(monA(:,3), monA(:,3))) monA(:,1) = monA(:,1)/sqrt(dot_product(monA(:,1), monA(:,1))) monA(:,2) = cross( monA(:,3), monA(:,1) ) write(*,'(/a/a)') . ' Local system of axes on A', . ' =========================' do i = 1, 3 write(*, '( 1x,i1,a,1x, 3(f12.7,2x) )' ) i,':', monA(:,i) enddo monB(:,3) = -(OH3 + OH4) monB(:,1) = (OH3 - OH4) monB(:,3) = monB(:,3)/sqrt(dot_product(monB(:,3), monB(:,3))) monB(:,1) = monB(:,1)/sqrt(dot_product(monB(:,1), monB(:,1))) monB(:,2) = cross( monB(:,3), monB(:,1) ) write(*,'(/a/a)') . ' Local system of axes on B', . ' =========================' do i = 1, 3 write(*, '( 1x,i1,a,1x, 3(f12.7,2x) )' ) i,':', monB(:,i) enddo Call Euler(monA, alphaA,betaA, gammaA) Call Euler(monB, alphaB,betaB, gammaB) alphaB = alphaB - alphaA alphaA = 0.d0 alphaA = alphaA*180.d0/pi betaA = betaA *180.d0/pi gammaA = gammaA*180.d0/pi alphaB = alphaB*180.d0/pi betaB = betaB *180.d0/pi gammaB = gammaB*180.d0/pi write(*,'(/a/a/ (a,f13.7) )') . ' Euler angles of A: ', . ' ================== ', . ' alpha = ', alphaA, . ' beta = ', betaA, . ' gamma = ', gammaA write(2,'(3 (f14.10,2x) )') alphaA, betaA, gammaA write(*,'(/a/a/ (a,f13.7) )') . ' Euler angles of B: ', . ' ================== ', . ' alpha = ', alphaB, . ' beta = ', betaB, . ' gamma = ', gammaB write(2,'(3 (f14.10,2x) )') alphaB, betaB, gammaB close(2) END Subroutine Euler(A, alpha,beta, gamma) C Determine Euler angles of the proper rotation matrix A. C Convention according to Biedenharn & Louck. use vector_functions, only: cart2polar implicit none intent(in) A intent(out) alpha, beta, gamma real*8 alpha, beta, gamma, . A(3,3), . polar3(3) C--Check if matrix is orthogonal with unit determinant. Call Check(A) C--Get polar angles of third column polar3 = cart2polar( A(:,3) ) ! returns R, theta, phi alpha = polar3(3) beta = polar3(2) C--Compute gamma Call Comgamma(A, alpha, beta, gamma) end Subroutine Check(A) C Check if matrix is orthogonal with unit determinant. use constants, only: mac_prec implicit none real*8 t, det, . A(3,3), . unit(3,3) intent(in) A unit = 0.d0 unit(1,1) = 1.d0 unit(2,2) = 1.d0 unit(3,3) = 1.d0 if (abs(sum(matmul(transpose(A),A)-unit)) .gt. mac_prec*1.d2 ) . then write(*,'(/a/a/3(d10.3) )') . ' Non-orthogonal matrix in subroutine check,', . ' overlap matrix: ', . matmul(transpose(A),A) endif t = det(A) if ( abs(t-1.d0) .gt. mac_prec*1.d2 ) then write(*,'(/a,1x,d14.6)') . ' Non-unit determinant:', t stop 16 endif end real*8 function det(A) C Determinant of A implicit none real*8 A(3,3), B(3) integer i intent(in) A B(1) = A(2,2)*A(3,3) - A(3,2)*A(2,3) B(2) = A(3,2)*A(1,3) - A(1,2)*A(3,3) B(3) = A(1,2)*A(2,3) - A(2,2)*A(1,3) det = 0.d0 do i=1,3 det = det + A(i,1)*B(i) enddo end Subroutine Comgamma(A, alpha, beta, gamma) C Compute third Euler angle gamma. use constants, only: pi, mac_prec implicit none intent(in) A, alpha, beta intent(out) gamma real*8 alpha, beta, gamma, . cb, sb, ca, sa, cg, sg, . A(3,3), . B1(3), . B2(3) integer i cb = cos(beta) sb = sin(beta) ca = cos(alpha) sa = sin(alpha) B1(1) = ca*cb B1(2) = sa*cb B1(3) = -sb B2(1) = -sa B2(2) = ca B2(3) = 0.d0 cg = 0.d0 sg = 0.d0 do i=1,3 cg = cg + B1(i)*A(i,1) sg = sg + B2(i)*A(i,1) enddo if ( cg .ge. 1.d0 ) then gamma = 0.d0 elseif ( cg .le. -1.d0 ) then gamma = pi else gamma = acos(cg) endif if ( sg < - mac_prec ) gamma = 2.d0*pi - gamma end