prg_syrotation_mod.F90

Go to the documentation of this file.

 
module prg_syrotation_mod
 
  implicit none
 
  private
 
  integer, parameter :: dp = kind(1.0d0)
 
  type, public :: rotation_type
    character(20) :: jobname
    character(50) :: typeofrot
    integer :: patom1
    integer :: patom2
    integer :: catom
    integer :: catom2
    real(dp) :: pq1(3)
    real(dp) :: pq2(3)
    real(dp) :: v1(3)
    real(dp) :: v2(3)
    real(dp) :: vq(3)
    integer :: rotate_atoms(2)
  end type rotation_type
 
  public :: prg_parse_rotation, prg_rotate
 
contains
 
  subroutine prg_parse_rotation(rot, filename)
 
    use prg_kernelparser_mod
    implicit none
    integer, parameter                ::  nkey_char = 2, nkey_int = 6, nkey_re = 15, nkey_log = 1
    character(20)                     ::  jobname
    character(len=*)                  ::  filename
    type(rotation_type), intent(inout)   ::  rot
 
    !Library of keywords with the respective defaults.
    character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
         'Jobname=','TypeOfRotation=']
    character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
         'MyJob','DirToDir']
 
    character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
         'Atom1=','Atom2=','CAtom=','CAtom2=','ListToRot1=','ListToRot2=']
    integer :: valvector_int(nkey_int) = (/ &
         0,0,0,0,1,10000000/)
 
    character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
         'PQ1X=','PQ1Y=','PQ1Z=','PQ2X=','PQ2Y=','PQ2Z=','V1X=','V1Y=', &
         &'V1Z=','V2X=','V2Y=','V2Z=','VQX=','VQY=','VQZ=']
    real(dp) :: valvector_re(nkey_re) = (/&
         0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/)
 
    character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
         'DUMMY=']
    logical :: valvector_log(nkey_log) = (/&
         .false./)
 
    !Start and stop characters
    character(len=50), parameter :: startstop(2) = [character(len=50) :: &
         'ROTATION{', '}']
 
    call prg_parsing_kernel(keyvector_char,valvector_char&
         ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
         keyvector_log,valvector_log,trim(filename),startstop)
 
    rot%jobname= valvector_char(1)
    rot%typeofrot= valvector_char(2)
 
    rot%patom1= valvector_int(1)
    rot%patom2= valvector_int(2)
    rot%catom= valvector_int(3)
    rot%catom2= valvector_int(4)
    rot%rotate_atoms(1)= valvector_int(5)
    rot%rotate_atoms(2)= valvector_int(6)
 
    rot%pq1(1)= valvector_re(1)
    rot%pq1(2)= valvector_re(2)
    rot%pq1(3)= valvector_re(3)
 
    rot%pq2(1)= valvector_re(4)
    rot%pq2(2)= valvector_re(5)
    rot%pq2(3)= valvector_re(6)
 
    rot%v1(1)= valvector_re(7)
    rot%v1(2)= valvector_re(8)
    rot%v1(3)= valvector_re(9)
 
    rot%v2(1)= valvector_re(10)
    rot%v2(2)= valvector_re(11)
    rot%v2(3)= valvector_re(12)
 
    rot%vQ(1)= valvector_re(13)
    rot%vQ(2)= valvector_re(14)
    rot%vQ(3)= valvector_re(15)
 
  end subroutine prg_parse_rotation
 
 
  subroutine prg_rotate(rot,r,verbose)
    integer                          ::  natoms, catom, catom2, i
    integer                          ::  k, l, patom1, patom2
    integer                          ::  rotate_atoms(2)
    integer, intent(in)              ::  verbose
    real(dp)                         ::  mat(3,3), p1minusc1(3), vtr(3)
    real(dp)                         ::  alpha, angle, d, deformation
    real(dp)                         ::  dotprod, dv1, dv2, hinge(3)
    real(dp)                         ::  hinge2, hingedotp1minusc1, pq1(3), pq2(3)
    real(dp)                         ::  pi, v1(3), v11(3), v1mod
    real(dp)                         ::  v1x, v1y, v1z, v1zversor
    real(dp)                         ::  v2(3), v22(3), vn(3), vq(3)
    real(dp)                         ::  vad(3), vdirn(3), vp(3), vr11(3)
    real(dp)                         ::  vr2(3), vr22(3)
    real(dp), allocatable            ::  rr(:,:), v(:,:), vr(:,:), vr1(:,:)
    real(dp), intent(inout)          ::  r(:,:)
    type(rotation_type), intent(in)  ::  rot
 
    pi=3.1415926535897993d0
    alpha= 0.0_dp
    alpha= (alpha/180.0_dp)*pi
 
    natoms = size(r, dim=2)
 
    allocate(v(3,natoms))
    allocate(vr(3,natoms))
    allocate(vr1(3,natoms))
    allocate(rr(3,natoms))
 
    vr=0.0_dp
    vr1=0.0_dp
    rr=0.0_dp
    v=0.0_dp
 
    patom1=rot%patom1           !Rotation atom indices
    patom2=rot%patom2
    catom=rot%catom       !Rotation center
    catom2=rot%catom2           !Second rotation center
 
    pq1(1)=rot%pq1(1)     !Initial director point
    pq1(2)=rot%pq1(2)
    pq1(3)=rot%pq1(3)
 
    pq2(1)=rot%pq2(1)     !Final director point
    pq2(2)=rot%pq2(2)
    pq2(3)=rot%pq2(3)
 
    v1(1)=rot%v1(1)         !Initial anchored vector
    v1(2)=rot%v1(2)
    v1(3)=rot%v1(3)
 
    v2(1)=rot%v2(1)         !Final anchored vector
    v2(2)=rot%v2(2)
    v2(3)=rot%v2(3)
 
    vq(1)=rot%vQ(1)           !Rotation center
    vq(2)=rot%vQ(2)
    vq(3)=rot%vQ(3)
 
    rotate_atoms(1)=rot%rotate_atoms(1)   !index for first rotated atom
    rotate_atoms(2)=rot%rotate_atoms(2)
 
    if(rotate_atoms(2).gt.natoms) rotate_atoms(2)=natoms
 
    if(patom1.ne.0)then
      pq1(1)=r(1,patom1)
      pq1(2)=r(2,patom1)
      pq1(3)=r(3,patom1)
 
      pq2(1)=r(1,patom2)
      pq2(2)=r(2,patom2)
      pq2(3)=r(3,patom2)
 
      vq(1)=r(1,catom)
      vq(2)=r(2,catom)
      vq(3)=r(3,catom)
    endif
 
    if(catom2.ne.0)then
      hinge=r(:,catom2)-r(:,catom)
      hinge2=hinge(1)**2+hinge(2)**2+hinge(3)**2
      p1minusc1=pq1(:)-r(:,catom)
      hingedotp1minusc1=hinge(1)*p1minusc1(1) + hinge(2)*p1minusc1(2) + hinge(3)*p1minusc1(3)
      vq(:)=r(:,catom) + hinge(:)*hingedotp1minusc1/hinge2
      write(*,*)vq(:)
      write(*,*)hinge(:)
      write(*,*)hingedotp1minusc1/hinge2
    endif
 
    if(v1(1).eq.0.0_dp.and.v1(2).eq.0.0_dp.and.v1(3).eq.0.0_dp)then
      v1=pq1-vq
    endif
 
    if(v2(1).eq.0.0_dp.and.v2(2).eq.0.0_dp.and.v2(3).eq.0.0_dp)then
      v2=pq2-vq
    endif
 
    vtr(1)=0.0_dp           !Translation
    vtr(2)=0.0_dp
    vtr(3)=0.0_dp
 
    dv1=dsqrt(v1(1)**2+v1(2)**2+v1(3)**2)  !Angle between v1 y v2
    dv2=dsqrt(v2(1)**2+v2(2)**2+v2(3)**2)
    dotprod=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
    alpha=-acos(dotprod/(dv1*dv2))
 
    if(alpha.eq.0.0_dp)write(*,*)'V1 and V2 are in the same direction'
    if(abs(alpha-pi).lt.1.0e-10)write(*,*)'V1 and V2 are in the same direction'
 
    if(verbose >= 1) write(*,*)'alpha=',alpha
 
    do i=rotate_atoms(1),rotate_atoms(2)
      r(1,i) = r(1,i)-vq(1)
      r(2,i) = r(2,i)-vq(2)
      r(3,i) = r(3,i)-vq(3)
    enddo
 
    vad=v1+v2
 
    vn(1)=v1(2)*v2(3)-v1(3)*v2(2)
    vn(2)=-v1(1)*v2(3)+v2(1)*v1(3)
    vn(3)=v1(1)*v2(2)-v2(1)*v1(2)
 
    vp(1)=vn(2)*vad(3)-vn(3)*vad(2)
    vp(2)=-vn(1)*vad(3)+vad(1)*vn(3)
    vp(3)=vn(1)*vad(2)-vad(1)*vn(2)
 
    d=dsqrt(vad(1)**2+vad(2)**2+vad(3)**2)
    vad=vad/d
 
    d=dsqrt(vp(1)**2+vp(2)**2+vp(3)**2)
    vp=vp/d
 
    d=dsqrt(vn(1)**2+vn(2)**2+vn(3)**2)
    vn=vn/d
 
    !Projection onto vad vp and vn
 
    do i=rotate_atoms(1),rotate_atoms(2)
      v(1,i)=r(1,i)*vad(1)+r(2,i)*vad(2)+r(3,i)*vad(3)
      v(2,i)=r(1,i)*vp(1)+r(2,i)*vp(2)+r(3,i)*vp(3)
      v(3,i)=r(1,i)*vn(1)+r(2,i)*vn(2)+r(3,i)*vn(3)
    enddo
 
    !Rotation matrix
 
    mat(1,1)=cos(alpha)
    mat(2,1)=sin(alpha)
    mat(3,1)=0.0_dp
    mat(1,2)=-sin(alpha)
    mat(2,2)=cos(alpha)
    mat(3,2)=0.0_dp
    mat(1,3)=0.0_dp
    mat(2,3)=0.0_dp
    mat(3,3)=1.0_dp
 
    !Rotation
 
    do l=rotate_atoms(1),rotate_atoms(2)
      do i=1,3
        do k=1,3
          vr(i,l)=vr(i,l)+mat(k,i)*v(k,l)
        enddo
      enddo
    enddo
 
    !Basis set change
 
    do i=rotate_atoms(1),rotate_atoms(2)
      vr1(1,i)=vr(1,i)*vad(1) + vr(2,i)*vp(1) + vr(3,i)*vn(1)
      vr1(2,i)=vr(1,i)*vad(2) + vr(2,i)*vp(2) + vr(3,i)*vn(2)
      vr1(3,i)=vr(1,i)*vad(3) + vr(2,i)*vp(3) + vr(3,i)*vn(3)
    enddo
 
    if(alpha.eq.0.0_dp)vr1=r
    if(abs(alpha-pi).lt.1.0e-10)vr1=-r
 
    do i=rotate_atoms(1),rotate_atoms(2)
      rr(1,i)=vr1(1,i)+vq(1)+vtr(1)
      rr(2,i)=vr1(2,i)+vq(2)+vtr(2)
      rr(3,i)=vr1(3,i)+vq(3)+vtr(3)
    enddo
 
    do i=rotate_atoms(1),rotate_atoms(2)
      r(1,i)=rr(1,i)
      r(2,i)=rr(2,i)
      r(3,i)=rr(3,i)
    enddo
 
    deformation=0.0d0
 
    do i=rotate_atoms(1),rotate_atoms(2)
      d=dsqrt((r(1,1)-r(1,i))**2+(r(2,1)-r(2,i))**2+(r(3,1)-r(3,i))**2)
      deformation=deformation + d
      d=dsqrt((rr(1,1)-rr(1,i))**2+(rr(2,1)-rr(2,i))**2+(rr(3,1)-rr(3,i))**2)
      deformation=deformation - d
    enddo
 
    if(deformation >= 0.001_dp)then
      write(*,*)'Rotation failed ...'
      write(*,*)'Deformation=',deformation
    endif
 
  end subroutine prg_rotate
 
end module prg_syrotation_mod