root / src / Tensor.f90 @ 4
Historique | Voir | Annoter | Télécharger (2,88 ko)
| 1 | 1 | equemene | subroutine tensor(natom,rx,ry,rz,Tinert,CoordG,rmass) |
|---|---|---|---|
| 2 | 1 | equemene | ! c |
| 3 | 1 | equemene | ! c elements of inertial tensor for a set of atomic systems. |
| 4 | 1 | equemene | ! c |
| 5 | 1 | equemene | ! c rx,ry,rz : arrays natom with cartesian coordinates of natom particles; input |
| 6 | 1 | equemene | ! c txx,tyy,tzz,txy,txz,tyz: elements of inertial tensors ; output |
| 7 | 1 | equemene | ! c a,b,c : returning center of mass coordinates. |
| 8 | 1 | equemene | ! c rmass : vector of length natom with particle masses; input |
| 9 | 1 | equemene | ! c natom : number of particles per system; input |
| 10 | 1 | equemene | ! c IOOUT : unit number for messages, silent for IOOUT.le.0; input |
| 11 | 1 | equemene | ! c subroutines called: cmshft m. lewerenz jul/92 |
| 12 | 1 | equemene | ! c |
| 13 | 1 | equemene | |
| 14 | 1 | equemene | use Io_module |
| 15 | 1 | equemene | |
| 16 | 1 | equemene | IMPLICIT NONE |
| 17 | 1 | equemene | |
| 18 | 1 | equemene | REAL(KREAL), parameter :: zero=0.d0,one=1.d0 |
| 19 | 1 | equemene | |
| 20 | 1 | equemene | INTEGER(KINT), INTENT(IN) :: Natom |
| 21 | 1 | equemene | real(KREAL) :: rx(Natom),ry(Natom),rz(Natom),rmass(Natom) |
| 22 | 1 | equemene | real(KREAL) :: a,b,c,Tinert(3,3),CoordG(3) |
| 23 | 1 | equemene | |
| 24 | 1 | equemene | INTEGER(KINT) :: jover, I, J |
| 25 | 1 | equemene | REAL(KREAL) :: Txx, Tyy, Tzz, Txy, Txz, Tyz, Txxi |
| 26 | 1 | equemene | REAL(KREAL) :: Ryb, Rzc, Rxa, Rmj2, ryb2, rzc2, rxa2, Rmj |
| 27 | 1 | equemene | |
| 28 | 1 | equemene | |
| 29 | 1 | equemene | if(natom.le.0) then |
| 30 | 1 | equemene | if(IOOUT.gt.0) write(IOOUT,'(/2a/)') & |
| 31 | 1 | equemene | ' *** error, illegal array dimension(s) in tensor,', & |
| 32 | 1 | equemene | ' return without action ***' |
| 33 | 1 | equemene | else if(natom.eq.1) then |
| 34 | 1 | equemene | a=rx(1) |
| 35 | 1 | equemene | b=ry(1) |
| 36 | 1 | equemene | c=rz(1) |
| 37 | 1 | equemene | txx=zero |
| 38 | 1 | equemene | tyy=zero |
| 39 | 1 | equemene | tzz=zero |
| 40 | 1 | equemene | txy=zero |
| 41 | 1 | equemene | txz=zero |
| 42 | 1 | equemene | tyz=zero |
| 43 | 1 | equemene | else |
| 44 | 1 | equemene | !c |
| 45 | 1 | equemene | !c first center of mass coordinates, then tensor elements |
| 46 | 1 | equemene | !c |
| 47 | 1 | equemene | call cmshft(natom,rx,ry,rz,rmass,a,b,c,0) |
| 48 | 1 | equemene | CoordG(1)=a |
| 49 | 1 | equemene | CoordG(2)=b |
| 50 | 1 | equemene | CoordG(3)=c |
| 51 | 1 | equemene | 1000 FORMAT(3F15.3) |
| 52 | 1 | equemene | !c WRITE(IOOUT,1000) (CoordG(i),i=1,3) |
| 53 | 1 | equemene | jover=mod(natom,2) |
| 54 | 1 | equemene | if(jover.eq.0) then |
| 55 | 1 | equemene | txx=zero |
| 56 | 1 | equemene | tyy=zero |
| 57 | 1 | equemene | tzz=zero |
| 58 | 1 | equemene | txy=zero |
| 59 | 1 | equemene | txz=zero |
| 60 | 1 | equemene | tyz=zero |
| 61 | 1 | equemene | else |
| 62 | 1 | equemene | rmj=rmass(1) |
| 63 | 1 | equemene | ryb=b-ry(1) |
| 64 | 1 | equemene | rzc=c-rz(1) |
| 65 | 1 | equemene | tyz=-rmj*ryb*rzc |
| 66 | 1 | equemene | rxa=a-rx(1) |
| 67 | 1 | equemene | txz=(-rmj*rxa)*rzc |
| 68 | 1 | equemene | txy=(-rmj*rxa)*ryb |
| 69 | 1 | equemene | txx=rmj*((ryb*ryb)+(rzc*rzc)) |
| 70 | 1 | equemene | tyy=rmj*((rxa*rxa)+(rzc*rzc)) |
| 71 | 1 | equemene | tzz=rmj*((rxa*rxa)+(ryb*ryb)) |
| 72 | 1 | equemene | end if |
| 73 | 1 | equemene | !c |
| 74 | 1 | equemene | !c force vectorization of inner loop if necessary (ibm-3090) |
| 75 | 1 | equemene | !c |
| 76 | 1 | equemene | do j=jover+1,natom,2 |
| 77 | 1 | equemene | rmj=rmass(j) |
| 78 | 1 | equemene | rmj2=rmass(j+1) |
| 79 | 1 | equemene | ryb=b-ry(j) |
| 80 | 1 | equemene | rzc=c-rz(j) |
| 81 | 1 | equemene | txxi=txx+rmj*((ryb*ryb)+(rzc*rzc)) |
| 82 | 1 | equemene | ryb2=b-ry(j+1) |
| 83 | 1 | equemene | rzc2=c-rz(j+1) |
| 84 | 1 | equemene | txx=txxi+rmj2*((ryb2*ryb2)+(rzc2*rzc2)) |
| 85 | 1 | equemene | tyz=tyz-rmj*ryb*rzc-rmj2*ryb2*rzc2 |
| 86 | 1 | equemene | rxa=a-rx(j) |
| 87 | 1 | equemene | rxa2=a-rx(j+1) |
| 88 | 1 | equemene | txy=txy-(rmj*rxa)*ryb-(rmj2*rxa2)*ryb2 |
| 89 | 1 | equemene | txz=txz-(rmj*rxa)*rzc-(rmj2*rxa2)*rzc2 |
| 90 | 1 | equemene | tyy=tyy+rmj*((rxa*rxa)+(rzc*rzc)) & |
| 91 | 1 | equemene | +rmj2*((rxa2*rxa2)+(rzc2*rzc2)) |
| 92 | 1 | equemene | tzz=tzz+rmj*((rxa*rxa)+(ryb*ryb)) & |
| 93 | 1 | equemene | +rmj2*((rxa2*rxa2)+(ryb2*ryb2)) |
| 94 | 1 | equemene | END DO |
| 95 | 1 | equemene | end if |
| 96 | 1 | equemene | Tinert(1,1)=txx |
| 97 | 1 | equemene | Tinert(1,2)=txy |
| 98 | 1 | equemene | Tinert(1,3)=txz |
| 99 | 1 | equemene | Tinert(2,1)=txy |
| 100 | 1 | equemene | Tinert(2,2)=tyy |
| 101 | 1 | equemene | Tinert(2,3)=tyz |
| 102 | 1 | equemene | Tinert(3,1)=txz |
| 103 | 1 | equemene | Tinert(3,2)=tyz |
| 104 | 1 | equemene | Tinert(3,3)=tzz |
| 105 | 1 | equemene | return |
| 106 | 1 | equemene | end subroutine tensor |