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