!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

