Chimie Théorique » Opt'n Path

!C================================================================

!C       Converti les positions Zmat en coordonnes cartesiennes

!C================================================================

        SUBROUTINE ConvertZmat_cart(iat,ind_zmat,val_zmat,x,y,z)

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

        use Path_module, only : Nat, KINT, KREAL

        IMPLICIT NONE

        integer(KINT) :: iat,n1,n2,n3

        real(KREAL)    :: x(nat),y(nat),z(nat)

        real(KREAL)    :: val_zmat(Nat,3)

        integer(KINT) :: ind_zmat(Nat,5)

        real(KREAL) ::  vx1,vy1,vz1,norm1

        real(KREAL) ::  vvx1,vvy1,vvz1,normv1

        real(KREAL) ::  vx2,vy2,vz2,norm2

        real(KREAL) ::  vvx2,vvy2,vvz2,normv2

        real(KREAL) ::  vx4, vy4, vz4

        real(KREAL) ::  d, a_val, a_dih

        real(KREAL) ::  a11_z1,a12_z1

        real(KREAL) ::  a11_z2,a12_z2

        real(KREAL) ::  a11_y,a12_y

!C ind_zmat(1) contient le numero de l'atome, (2) celui par rapport auquel on definit la distance...

          n1=ind_zmat(iat,2)

          n2=ind_zmat(iat,3)

          n3=ind_zmat(iat,4)

!          WRITE(*,*) iat,n1,n2,n3

          d=val_zmat(iat,1)

          a_val=val_zmat(iat,2)/180.*3.141592654

          a_dih=val_zmat(iat,3)/180.*3.141592654

!          WRITE(*,*) "d,val,di",d,a_val,a_dih,a_val*180./3.141592654

          CALL vecteur(n1,n2,x,y,z,vx1,vy1,vz1,norm1)

          CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)

          vvx1=vx1

          vvy1=vy1

          vvz1=vz1

!c rotation autour de z de v1 et v2 de phi (a11_z1=cos(phi) et a12_z1=sin(phi) ) tq

!C  v1 soit dans le plan ????

          normv1=dsqrt(vvx1*vvx1+vvy1*vvy1)

          IF (normv1 .GE. 1.D-6) THEN

             a11_z1 = vvx1/normv1

             a12_z1 = vvy1/normv1

           ELSE

             a11_z1 = 1

             a12_z1 = 0

          END IF

             CALL rota_z(vx1,vy1,vz1,a11_z1,-a12_z1)

             CALL rota_z(vx2,vy2,vz2,a11_z1,-a12_z1)

!c rotation autour de y de v1 et v2 de theta (a11_y=cos(theta) et a12_y=sin(theta) ) tq

!C  v1 soit dans le plan ????

          IF (norm1 .GE. 1.D-8) THEN

             a11_y = vz1/norm1

             a12_y = vx1/norm1

          ELSE

             a11_y = 1

             a12_y = 0

          END IF

             CALL rota_y(vx1,vy1,vz1,a11_y,a12_y)

             CALL rota_y(vx2,vy2,vz2,a11_y,a12_y)

!c rotation autour de z de v1 et v2 de psi (a11_z2=cos(psi) et a12_z2=sin(psi) ) tq

!C  v1 soit dans le plan ????

          vvx2=vx2

          vvy2=vy2

          vvz2=vz2

          normv2=dsqrt(vvx2*vvx2+vvy2*vvy2)

          IF (normv2 .GE. 1.D-8) THEN

             a11_z2 = vvx2/normv2

             a12_z2 = vvy2/normv2

          ELSE

             a11_z2 = 1

             a12_z2 = 0

          END IF

             CALL rota_z(vx1,vy1,vz1,a11_z2,-a12_z2)

             CALL rota_z(vx2,vy2,vz2,a11_z2,-a12_z2)

!c calcul le vecteur de l atome dans la nouvelle orientation

          vx4=d*dsin(a_val)*dcos(a_dih)

          vy4=-d*dsin(a_val)*dsin(a_dih)

           vz4=d*cos(a_val)

!c calcul le vecteur de l atome

!c on tourne en sens invers

          CALL rota_z(vx4,vy4,vz4,a11_z2, a12_z2)

          CALL rota_y(vx4,vy4,vz4,a11_y ,-a12_y )

          CALL rota_z(vx4,vy4,vz4,a11_z1, a12_z1)

!c calcul l atome a partir de v4

          x(iat)=vx4+x(n1)

          y(iat)=vy4+y(n1)

          z(iat)=vz4+z(n1)

        END