root / src / ConvertZmat_cart.f90 @ 1
Historique | Voir | Annoter | Télécharger (3,32 ko)
| 1 | 1 | equemene | !C================================================================ |
|---|---|---|---|
| 2 | 1 | equemene | !C Converti les positions Zmat en coordonnes cartesiennes |
| 3 | 1 | equemene | !C================================================================ |
| 4 | 1 | equemene | |
| 5 | 1 | equemene | SUBROUTINE ConvertZmat_cart(iat,ind_zmat,val_zmat,x,y,z) |
| 6 | 1 | equemene | |
| 7 | 1 | equemene | use Path_module, only : Nat, KINT, KREAL |
| 8 | 1 | equemene | |
| 9 | 1 | equemene | IMPLICIT NONE |
| 10 | 1 | equemene | |
| 11 | 1 | equemene | integer(KINT) :: iat,n1,n2,n3 |
| 12 | 1 | equemene | real(KREAL) :: x(nat),y(nat),z(nat) |
| 13 | 1 | equemene | real(KREAL) :: val_zmat(Nat,3) |
| 14 | 1 | equemene | integer(KINT) :: ind_zmat(Nat,5) |
| 15 | 1 | equemene | |
| 16 | 1 | equemene | real(KREAL) :: vx1,vy1,vz1,norm1 |
| 17 | 1 | equemene | real(KREAL) :: vvx1,vvy1,vvz1,normv1 |
| 18 | 1 | equemene | real(KREAL) :: vx2,vy2,vz2,norm2 |
| 19 | 1 | equemene | real(KREAL) :: vvx2,vvy2,vvz2,normv2 |
| 20 | 1 | equemene | real(KREAL) :: vx3,vy3,vz3,norm3 |
| 21 | 1 | equemene | real(KREAL) :: vx4,vy4,vz4,norm4 |
| 22 | 1 | equemene | real(KREAL) :: d,a_val,a_dih,c1,c2 |
| 23 | 1 | equemene | real(KREAL) :: a11_z1,a12_z1 |
| 24 | 1 | equemene | real(KREAL) :: a11_z2,a12_z2 |
| 25 | 1 | equemene | real(KREAL) :: a11_y,a12_y |
| 26 | 1 | equemene | |
| 27 | 1 | equemene | |
| 28 | 1 | equemene | |
| 29 | 1 | equemene | !C ind_zmat(1) contient le numero de l'atome, (2) celui par rapport auquel on definit la distance... |
| 30 | 1 | equemene | |
| 31 | 1 | equemene | n1=ind_zmat(iat,2) |
| 32 | 1 | equemene | n2=ind_zmat(iat,3) |
| 33 | 1 | equemene | n3=ind_zmat(iat,4) |
| 34 | 1 | equemene | |
| 35 | 1 | equemene | ! WRITE(*,*) iat,n1,n2,n3 |
| 36 | 1 | equemene | |
| 37 | 1 | equemene | d=val_zmat(iat,1) |
| 38 | 1 | equemene | a_val=val_zmat(iat,2)/180.*3.141592654 |
| 39 | 1 | equemene | a_dih=val_zmat(iat,3)/180.*3.141592654 |
| 40 | 1 | equemene | |
| 41 | 1 | equemene | ! WRITE(*,*) "d,val,di",d,a_val,a_dih,a_val*180./3.141592654 |
| 42 | 1 | equemene | |
| 43 | 1 | equemene | CALL vecteur(n1,n2,x,y,z,vx1,vy1,vz1,norm1) |
| 44 | 1 | equemene | CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2) |
| 45 | 1 | equemene | |
| 46 | 1 | equemene | vvx1=vx1 |
| 47 | 1 | equemene | vvy1=vy1 |
| 48 | 1 | equemene | vvz1=vz1 |
| 49 | 1 | equemene | |
| 50 | 1 | equemene | !c rotation autour de z de v1 et v2 de phi (a11_z1=cos(phi) et a12_z1=sin(phi) ) tq |
| 51 | 1 | equemene | !C v1 soit dans le plan ???? |
| 52 | 1 | equemene | |
| 53 | 1 | equemene | normv1=dsqrt(vvx1*vvx1+vvy1*vvy1) |
| 54 | 1 | equemene | |
| 55 | 1 | equemene | IF (normv1 .GE. 1.D-6) THEN |
| 56 | 1 | equemene | a11_z1 = vvx1/normv1 |
| 57 | 1 | equemene | a12_z1 = vvy1/normv1 |
| 58 | 1 | equemene | ELSE |
| 59 | 1 | equemene | a11_z1 = 1 |
| 60 | 1 | equemene | a12_z1 = 0 |
| 61 | 1 | equemene | END IF |
| 62 | 1 | equemene | |
| 63 | 1 | equemene | CALL rota_z(vx1,vy1,vz1,a11_z1,-a12_z1) |
| 64 | 1 | equemene | CALL rota_z(vx2,vy2,vz2,a11_z1,-a12_z1) |
| 65 | 1 | equemene | |
| 66 | 1 | equemene | |
| 67 | 1 | equemene | !c rotation autour de y de v1 et v2 de theta (a11_y=cos(theta) et a12_y=sin(theta) ) tq |
| 68 | 1 | equemene | !C v1 soit dans le plan ???? |
| 69 | 1 | equemene | |
| 70 | 1 | equemene | |
| 71 | 1 | equemene | IF (norm1 .GE. 1.D-8) THEN |
| 72 | 1 | equemene | a11_y = vz1/norm1 |
| 73 | 1 | equemene | a12_y = vx1/norm1 |
| 74 | 1 | equemene | ELSE |
| 75 | 1 | equemene | a11_y = 1 |
| 76 | 1 | equemene | a12_y = 0 |
| 77 | 1 | equemene | END IF |
| 78 | 1 | equemene | |
| 79 | 1 | equemene | CALL rota_y(vx1,vy1,vz1,a11_y,a12_y) |
| 80 | 1 | equemene | CALL rota_y(vx2,vy2,vz2,a11_y,a12_y) |
| 81 | 1 | equemene | |
| 82 | 1 | equemene | |
| 83 | 1 | equemene | !c rotation autour de z de v1 et v2 de psi (a11_z2=cos(psi) et a12_z2=sin(psi) ) tq |
| 84 | 1 | equemene | !C v1 soit dans le plan ???? |
| 85 | 1 | equemene | |
| 86 | 1 | equemene | vvx2=vx2 |
| 87 | 1 | equemene | vvy2=vy2 |
| 88 | 1 | equemene | vvz2=vz2 |
| 89 | 1 | equemene | normv2=dsqrt(vvx2*vvx2+vvy2*vvy2) |
| 90 | 1 | equemene | |
| 91 | 1 | equemene | IF (normv2 .GE. 1.D-8) THEN |
| 92 | 1 | equemene | a11_z2 = vvx2/normv2 |
| 93 | 1 | equemene | a12_z2 = vvy2/normv2 |
| 94 | 1 | equemene | ELSE |
| 95 | 1 | equemene | a11_z2 = 1 |
| 96 | 1 | equemene | a12_z2 = 0 |
| 97 | 1 | equemene | END IF |
| 98 | 1 | equemene | |
| 99 | 1 | equemene | CALL rota_z(vx1,vy1,vz1,a11_z2,-a12_z2) |
| 100 | 1 | equemene | CALL rota_z(vx2,vy2,vz2,a11_z2,-a12_z2) |
| 101 | 1 | equemene | |
| 102 | 1 | equemene | |
| 103 | 1 | equemene | !c calcul le vecteur de l atome dans la nouvelle orientation |
| 104 | 1 | equemene | |
| 105 | 1 | equemene | vx4=d*dsin(a_val)*dcos(a_dih) |
| 106 | 1 | equemene | vy4=-d*dsin(a_val)*dsin(a_dih) |
| 107 | 1 | equemene | vz4=d*cos(a_val) |
| 108 | 1 | equemene | |
| 109 | 1 | equemene | !c calcul le vecteur de l atome |
| 110 | 1 | equemene | !c on tourne en sens invers |
| 111 | 1 | equemene | CALL rota_z(vx4,vy4,vz4,a11_z2, a12_z2) |
| 112 | 1 | equemene | CALL rota_y(vx4,vy4,vz4,a11_y ,-a12_y ) |
| 113 | 1 | equemene | CALL rota_z(vx4,vy4,vz4,a11_z1, a12_z1) |
| 114 | 1 | equemene | |
| 115 | 1 | equemene | !c calcul l atome a partir de v4 |
| 116 | 1 | equemene | x(iat)=vx4+x(n1) |
| 117 | 1 | equemene | y(iat)=vy4+y(n1) |
| 118 | 1 | equemene | z(iat)=vz4+z(n1) |
| 119 | 1 | equemene | |
| 120 | 1 | equemene | END |