Révision 12
| src/Std_ori.f90 (revision 12) | ||
|---|---|---|
| 1 |
subroutine Std_ori(natom,rx,ry,rz,a,b,c,rmass) |
|
| 2 |
!c |
|
| 3 |
!c rotational constants for all replicas by inline diagonalization |
|
| 4 |
!c of inertial tensor |
|
| 5 |
!c |
|
| 6 |
!!!!!!! Input |
|
| 7 |
! natom : number of atoms |
|
| 8 |
! rx,ry,rz : arrays natom with cartesian coordinates of natom particles |
|
| 9 |
! rmass : vector of length natom with particle masses |
|
| 10 |
! |
|
| 11 |
!!!!!!! Output |
|
| 12 |
! a,b,c : rotational constants a, b, and c |
|
| 13 |
! |
|
| 14 |
!!!!!! Dependencies |
|
| 15 |
! tensor: calculate the inertial tensor |
|
| 16 |
! Jacobi, trie : diaganalization subroutines |
|
| 17 |
|
|
| 18 |
implicit none |
|
| 19 |
|
|
| 20 |
INTEGER, PARAMETER :: KINT=KIND(1) |
|
| 21 |
INTEGER, PARAMETER :: KREAL=KIND(1.0D0) |
|
| 22 |
|
|
| 23 |
INTEGER(KINT), INTENT(IN) :: Natom |
|
| 24 |
REAL(KREAL), INTENT(INOUT) :: rx(Natom),ry(Natom),rz(Natom) |
|
| 25 |
REAL(KREAL), INTENT(IN) :: rmass(Natom) |
|
| 26 |
REAL(KREAL), INTENT(OUT) :: a,b,c |
|
| 27 |
|
|
| 28 |
real(KREAL) :: Tinert(3,3),D(3),V(3,3),CoordG(3) |
|
| 29 |
real(KREAL) :: ca,sa,na,cb,sb,nb |
|
| 30 |
|
|
| 31 |
INTEGER(KINT) :: I |
|
| 32 |
REAL(KREAL) :: Vt,Rxt,Ryt |
|
| 33 |
|
|
| 34 |
call tensor(natom,rx,ry,rz,Tinert,CoordG,rmass) |
|
| 35 |
call Jacobi(Tinert,3,D,V,3) |
|
| 36 |
call Trie(3,D,V,3) |
|
| 37 |
|
|
| 38 |
!1000 FORMAT(3F15.3) |
|
| 39 |
!1010 FORMAT(4F15.3) |
|
| 40 |
|
|
| 41 |
a=D(1) |
|
| 42 |
! WRITE(*,1010) a,(V(i,1),i=1,3) |
|
| 43 |
b=D(2) |
|
| 44 |
! WRITE(*,1010) b,(V(i,2),i=1,3) |
|
| 45 |
c=D(3) |
|
| 46 |
! WRITE(*,1010) c,(V(i,3),i=1,3) |
|
| 47 |
|
|
| 48 |
!C On commence par mettre le centre de masse a l'origine |
|
| 49 |
|
|
| 50 |
Do I=1,natom |
|
| 51 |
! WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I) |
|
| 52 |
rx(i)=rx(i)-CoordG(1) |
|
| 53 |
ry(i)=ry(i)-CoordG(2) |
|
| 54 |
rz(i)=rz(i)-CoordG(3) |
|
| 55 |
! WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I) |
|
| 56 |
END DO |
|
| 57 |
|
|
| 58 |
! C Ensuite, on va faire coincider les axes d'inertie avec x,y,z |
|
| 59 |
! C Pour etre coherent, et bien que ce ne soit pas le plus logique, |
|
| 60 |
! C on classe les vecteurs dans le sens croissant, et on fait coincider |
|
| 61 |
! C les axes dans l'ordre x,y,z pour l'instant. |
|
| 62 |
|
|
| 63 |
! C Coincidation du premier axe avec x |
|
| 64 |
! C Rotation autour de z pour l'amener dans le plan xOz |
|
| 65 |
! C puis Rotation autour de y pour l'amener suivant x |
|
| 66 |
na=sqrt(V(1,1)*V(1,1)+V(2,1)*V(2,1)) |
|
| 67 |
nb=1 |
|
| 68 |
if (na.gt.1D-10) then |
|
| 69 |
ca=V(1,1)/na |
|
| 70 |
sa=-V(2,1)/na |
|
| 71 |
else |
|
| 72 |
ca=1 |
|
| 73 |
sa=0 |
|
| 74 |
END IF |
|
| 75 |
cb=na/nb |
|
| 76 |
sb=-V(3,1)/nb |
|
| 77 |
DO I=1,natom |
|
| 78 |
rxt=rx(i) |
|
| 79 |
rx(i)=rx(i)*ca-ry(i)*sa |
|
| 80 |
ry(i)=rxt*sa+ry(i)*ca |
|
| 81 |
rxt=rx(i) |
|
| 82 |
rx(i)=rx(i)*cb-rz(i)*sb |
|
| 83 |
rz(i)=rxt*sb+rz(i)*cb |
|
| 84 |
END DO |
|
| 85 |
|
|
| 86 |
!C on applique aussi la rotation aux deux autres axes... |
|
| 87 |
Vt=V(1,2) |
|
| 88 |
V(1,2)=V(1,2)*ca-V(2,2)*sa |
|
| 89 |
V(2,2)=Vt*sa+V(2,2)*ca |
|
| 90 |
Vt=V(1,2) |
|
| 91 |
V(1,2)=V(1,2)*cb-V(3,2)*sb |
|
| 92 |
V(3,2)=Vt*sb+V(3,2)*cb |
|
| 93 |
|
|
| 94 |
!c WRITE(IOOUT,1010) b,(V(i,2),i=1,3) |
|
| 95 |
|
|
| 96 |
!c Vt=V(1,3) |
|
| 97 |
!c V(1,3)=V(1,3)*ca-V(2,3)*sa |
|
| 98 |
!c V(2,3)=Vt*sa+V(2,3)*ca |
|
| 99 |
!c Vt=V(1,3) |
|
| 100 |
!c V(1,3)=V(1,3)*cb-V(3,3)*sb |
|
| 101 |
!c V(3,3)=Vt*sb+V(3,3)*cb |
|
| 102 |
!c WRITE(IOOUT,1010) c,(V(i,3),i=1,3) |
|
| 103 |
|
|
| 104 |
!C et on repart pour le 2e axe sur y |
|
| 105 |
!C commes les vecteurs sont orthonormes, il est |
|
| 106 |
!C deja dans le plan yOz d'ou une seule rotation |
|
| 107 |
!c na=sqrt(V(2,2)*V(2,2)+V(3,2)*V(3,2)) |
|
| 108 |
na=1 |
|
| 109 |
nb=1 |
|
| 110 |
ca=V(2,2)/na |
|
| 111 |
sa=-V(3,2)/na |
|
| 112 |
!c write(IOOUT,*) 'ca:',ca,sa,cb,sb,na |
|
| 113 |
DO I=1,natom |
|
| 114 |
ryt=ry(i) |
|
| 115 |
ry(i)=ry(i)*ca-rz(i)*sa |
|
| 116 |
rz(i)=ryt*sa+rz(i)*ca |
|
| 117 |
END DO |
|
| 118 |
!C on applique aussi la rotation au dernier axe... |
|
| 119 |
!c Vt=V(2,3) |
|
| 120 |
!c V(2,3)=V(2,3)*ca-V(3,3)*sa |
|
| 121 |
!c V(3,3)=Vt*sa+V(3,3)*ca |
|
| 122 |
!c WRITE(IOOUT,1010) c,(V(i,3),i=1,3) |
|
| 123 |
|
|
| 124 |
!C et le 3e axe est deja sur z car ils sont orthogonaux... |
|
| 125 |
|
|
| 126 |
return |
|
| 127 |
end subroutine Std_ori |
|
| src/egrad_chamfre.f90 (revision 12) | ||
|---|---|---|
| 1 |
!*************############### Attention ###############************* |
|
| 2 |
!(a) Two parameters aa and br(3,3) must be provided in main program |
|
| 3 |
!(b) At first, z=0 for the outmost surface |
|
| 4 |
!(c) The units of all parameters are not atomic units, but angstrom, eV etc. |
|
| 5 |
!(d) Before using this routine, you must confirm that the energy is CONSISTENT with its derivative. |
|
| 6 |
! This check of consistency is VERY IMPORTANT for energy conservation in MD. |
|
| 7 |
!!! |
|
| 8 |
!! xa, ya,za present the positions of atoms. v is the energy of the system. dvdr present the force along the x,y,z directions. |
|
| 9 |
!! |
|
| 10 |
subroutine egrad_chamfre(Nat,V,Geom,Grad) |
|
| 11 |
! WAS xa,ya,za,v,dvdr) |
|
| 12 |
|
|
| 13 |
|
|
| 14 |
use Path_module, only : Lat_a,Lat_b,order,renum |
|
| 15 |
use Io_module |
|
| 16 |
|
|
| 17 |
|
|
| 18 |
IMPLICIT NONE |
|
| 19 |
|
|
| 20 |
INTEGER(KINT), INTENT(IN) :: Nat |
|
| 21 |
REAL(KREAL), INTENT(OUT) :: V,grad(Nat*3) |
|
| 22 |
REAL(KREAL), INTENT(IN) :: Geom(Nat,3) |
|
| 23 |
|
|
| 24 |
! Bohr --> Angstr |
|
| 25 |
! real(KREAL), parameter :: BOHR = 0.52917726D+00 |
|
| 26 |
|
|
| 27 |
!! PFL 2010 Nov. 28 -> |
|
| 28 |
! for now I keep this old F77 scheme. This has to be replaced by ALLOCATABLE |
|
| 29 |
! things |
|
| 30 |
! <- PFL 2010 Nov 28. |
|
| 31 |
REAL(KREAL) :: xa(100),ya(100),za(100),dvdr(100,3) |
|
| 32 |
REAL(KREAL) :: rshort,rlong,zshort,zlong |
|
| 33 |
REAL(KREAL), SAVE :: rnc, rnn, rhh, rch, rnh, ann(100), anc(100), anh(100), ach(100), ahh(100) |
|
| 34 |
REAL(KREAL) :: dv(100,10,3),df(100,10,3) |
|
| 35 |
REAL(KREAL) :: px(9000),py(9000),pz(9000),tmp,temp,rtemp,yf |
|
| 36 |
INTEGER(KINT) :: ka(9000),kx(9000),ky(9000),itx,ity,nf(100),mf(100,500) |
|
| 37 |
REAL(KREAL) :: xi(3), xj(3), rij, er, eb |
|
| 38 |
REAL(KREAL) :: fcut, tf |
|
| 39 |
INTEGER(KINT) :: np,nq,nh,nc |
|
| 40 |
INTEGER(KINT) :: i, j, nn, m, n, mid, n1, im, jm, n2,n3,n4, iat |
|
| 41 |
LOGICAL, SAVE :: First=.TRUE. |
|
| 42 |
LOGICAL :: debug, FExist |
|
| 43 |
|
|
| 44 |
common/zdis/zshort,zlong |
|
| 45 |
common/rdistanc/rshort,rlong |
|
| 46 |
|
|
| 47 |
|
|
| 48 |
INTERFACE |
|
| 49 |
function valid(string) result (isValid) |
|
| 50 |
CHARACTER(*), intent(in) :: string |
|
| 51 |
logical :: isValid |
|
| 52 |
END function VALID |
|
| 53 |
|
|
| 54 |
END INTERFACE |
|
| 55 |
|
|
| 56 |
!********************************************** |
|
| 57 |
|
|
| 58 |
|
|
| 59 |
debug=valid('egrad_chamfre')
|
|
| 60 |
|
|
| 61 |
if (debug) Call Header('Entering Egrad_chamfre')
|
|
| 62 |
|
|
| 63 |
v=0.0;dvdr=0.0 |
|
| 64 |
!******prepare the parameters of REBO potential |
|
| 65 |
!rshort=3.5 |
|
| 66 |
!rlong=3.9 |
|
| 67 |
zshort=3.5 |
|
| 68 |
zlong=4.5 |
|
| 69 |
|
|
| 70 |
rnn=3.2 !! to adjust |
|
| 71 |
rnc=3.5 !! to adjust |
|
| 72 |
rnh=3.0 !! to adjust |
|
| 73 |
rch=3.2 !! to adjust |
|
| 74 |
rhh=3.2 !! to adjust |
|
| 75 |
|
|
| 76 |
!rph=rlong |
|
| 77 |
|
|
| 78 |
if (First) THEN |
|
| 79 |
First=.FALSE. |
|
| 80 |
INQUIRE(FILE="parameter.dat",EXIST=FExist) |
|
| 81 |
if (.NOT.Fexist) THEN |
|
| 82 |
WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!" |
|
| 83 |
WRITE(*,*) "!! egrad_chamfre !!!!" |
|
| 84 |
WRITE(*,*) "The file parameter.dat does not exists." |
|
| 85 |
WRITE(*,*) "!! egrad_chamfre !!!!" |
|
| 86 |
WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!" |
|
| 87 |
STOP |
|
| 88 |
END IF |
|
| 89 |
open(1,file="parameter.dat") |
|
| 90 |
do i=1,5 !Ni-Ni |
|
| 91 |
read(1,*)ann(i) |
|
| 92 |
enddo |
|
| 93 |
|
|
| 94 |
do i=1,5 !Ni-C |
|
| 95 |
read(1,*)anc(i) |
|
| 96 |
enddo |
|
| 97 |
|
|
| 98 |
do i=1,5 !Ni-H |
|
| 99 |
read(1,*)anh(i) |
|
| 100 |
enddo |
|
| 101 |
|
|
| 102 |
do i=1,5 !H-H |
|
| 103 |
read(1,*)ahh(i) |
|
| 104 |
enddo |
|
| 105 |
do i=1,5 !C-H |
|
| 106 |
read(1,*)ach(i) |
|
| 107 |
enddo |
|
| 108 |
!********LONG INTERACTION ********** |
|
| 109 |
!do i=1,2 |
|
| 110 |
!read(1,*)a(i) |
|
| 111 |
!enddo |
|
| 112 |
close(1) |
|
| 113 |
END IF |
|
| 114 |
!************************************ |
|
| 115 |
np=45 |
|
| 116 |
nh=4 |
|
| 117 |
nc=1 |
|
| 118 |
nq=np+nh+nc |
|
| 119 |
!************************************ |
|
| 120 |
|
|
| 121 |
DO I=1,Nat |
|
| 122 |
Iat=Order(I) |
|
| 123 |
xa(I)=Geom(Iat,1) |
|
| 124 |
ya(I)=Geom(Iat,2) |
|
| 125 |
za(I)=Geom(Iat,3) |
|
| 126 |
END DO |
|
| 127 |
|
|
| 128 |
px=0.0;py=0.0;pz=0.0 |
|
| 129 |
nn=0 |
|
| 130 |
do m=-1,1 |
|
| 131 |
do n=-1,1 |
|
| 132 |
do i=1,nq |
|
| 133 |
if(m==0.and.n==0.and.i==1)mid=nn |
|
| 134 |
nn=nn+1 |
|
| 135 |
px(nn)=xa(i)+m*lat_a(1)+n*lat_b(1) |
|
| 136 |
py(nn)=ya(i)+m*lat_a(2)+n*lat_b(2) |
|
| 137 |
pz(nn)=za(i) |
|
| 138 |
ka(nn)=i;kx(nn)=m;ky(nn)=n |
|
| 139 |
enddo |
|
| 140 |
enddo |
|
| 141 |
enddo |
|
| 142 |
!***********************xyz file required for MS input |
|
| 143 |
open(123,file="111.xyz",ACCESS='APPEND') |
|
| 144 |
write(123,*)nn |
|
| 145 |
write(123,*) |
|
| 146 |
do i=1,nn |
|
| 147 |
if(ka(i)<(np+1))then |
|
| 148 |
write(123,10)'Ni',px(i),py(i),pz(i) |
|
| 149 |
elseif (ka(i)<np+5) then |
|
| 150 |
write(123,10)'H',px(i),py(i),pz(i) |
|
| 151 |
else |
|
| 152 |
write(123,10)'C',px(i),py(i),pz(i) |
|
| 153 |
10 format(a6,3f15.8) |
|
| 154 |
endif |
|
| 155 |
enddo |
|
| 156 |
close(123) |
|
| 157 |
!****************************************** |
|
| 158 |
do i=1,nq |
|
| 159 |
tmp=px(i);px(i)=px(i+mid);px(i+mid)=tmp |
|
| 160 |
tmp=py(i);py(i)=py(i+mid);py(i+mid)=tmp |
|
| 161 |
tmp=pz(i);pz(i)=pz(i+mid);pz(i+mid)=tmp |
|
| 162 |
itx=kx(i);kx(i)=kx(i+mid);kx(i+mid)=itx |
|
| 163 |
ity=ky(i);ky(i)=ky(i+mid);ky(i+mid)=ity |
|
| 164 |
enddo |
|
| 165 |
mid=0 |
|
| 166 |
!****************************************** |
|
| 167 |
nf=0;mf=0 |
|
| 168 |
do i=1,nq |
|
| 169 |
nf(i)=0 |
|
| 170 |
do j=1,nn |
|
| 171 |
temp=sqrt((px(i)-px(j))**2+(py(i)-py(j))**2+(pz(i)-pz(j))**2) |
|
| 172 |
! |
|
| 173 |
!rtemp=rlong |
|
| 174 |
if (ka(i+mid)<np+1.and.ka(j+mid)<np+1)then ! Ni-Ni |
|
| 175 |
|
|
| 176 |
rtemp=rnn |
|
| 177 |
endif |
|
| 178 |
if (ka(i+mid)<np+1.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then ! Ni-H |
|
| 179 |
|
|
| 180 |
rtemp=rnh |
|
| 181 |
endif |
|
| 182 |
if (ka(i+mid)<np+1.and.ka(j+mid)>np+4) then ! Ni-C |
|
| 183 |
|
|
| 184 |
rtemp=rnc |
|
| 185 |
endif |
|
| 186 |
if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np+4) then ! H-C |
|
| 187 |
|
|
| 188 |
rtemp=rch |
|
| 189 |
endif |
|
| 190 |
if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then ! H-H |
|
| 191 |
|
|
| 192 |
rtemp=rhh |
|
| 193 |
endif |
|
| 194 |
if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)<np+1) then ! H-Ni |
|
| 195 |
|
|
| 196 |
rtemp=rnh |
|
| 197 |
endif |
|
| 198 |
if (ka(i+mid)>np+4.and.ka(j+mid)<np+1) then ! C-Ni |
|
| 199 |
|
|
| 200 |
rtemp=rnc |
|
| 201 |
endif |
|
| 202 |
if (ka(i+mid)>np+4.and.ka(j+mid)>np.and.ka(j+mid)<np+5) then ! C-H |
|
| 203 |
|
|
| 204 |
rtemp=rch |
|
| 205 |
endif |
|
| 206 |
|
|
| 207 |
if(temp>0.1.and.temp<rtemp)then |
|
| 208 |
nf(i)=nf(i)+1 |
|
| 209 |
mf(i,nf(i))=j |
|
| 210 |
endif |
|
| 211 |
! |
|
| 212 |
enddo |
|
| 213 |
!write(*,*)nf(i) |
|
| 214 |
enddo |
|
| 215 |
!************************************** |
|
| 216 |
!stop |
|
| 217 |
|
|
| 218 |
!*************************************Starting to calculate the sum |
|
| 219 |
yf=0.0;dvdr=0.0 |
|
| 220 |
do i=1,nq |
|
| 221 |
er=0.0;eb=0.0 |
|
| 222 |
dv=0.0;df=0.0 |
|
| 223 |
do j=1,nf(i) |
|
| 224 |
|
|
| 225 |
xi(1)=px(i); xi(2)=py(i); xi(3)=pz(i) |
|
| 226 |
xj(1)=px(mf(i,j));xj(2)=py(mf(i,j));xj(3)=pz(mf(i,j)) |
|
| 227 |
rij=sqrt((xj(1)-xi(1))**2.0+(xj(2)-xi(2))**2.0+(xj(3)-xi(3))**2.0) |
|
| 228 |
n1=ka(i)/(np+1);n2=ka(mf(i,j))/(np+1) !!! Ni 46 H 50 C |
|
| 229 |
n3=ka(i)/(np+5);n4=ka(mf(i,j))/(np+5) |
|
| 230 |
|
|
| 231 |
!**********************************Ni-Ni interaction |
|
| 232 |
! ********** |
|
| 233 |
if(n1==0.and.n3==0.and.n2==0.and.n4==0)then !Ni-Ni |
|
| 234 |
call fctNiNi(rnn,rij,fcut) |
|
| 235 |
call dfcutNiNi(rnn,rij,tf) |
|
| 236 |
eb=eb+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut |
|
| 237 |
er=er+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut |
|
| 238 |
|
|
| 239 |
!Calculation of FORCE |
|
| 240 |
do im=1,3 |
|
| 241 |
dv(i,1,im)=dv(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut & |
|
| 242 |
*(-ann(2)/ann(3))*(xi(im)-xj(im))/rij |
|
| 243 |
dv(i,2,im)=dv(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut |
|
| 244 |
dv(i,3,im)=dv(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut & |
|
| 245 |
*(-2.0*ann(1)/ann(3))*(xi(im)-xj(im))/rij |
|
| 246 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) & |
|
| 247 |
*(rij/ann(3)-1.0))*fcut*(-ann(2)/ann(3))*(xj(im)-xi(im))/rij |
|
| 248 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) & |
|
| 249 |
*(rij/ann(3)-1.0))*fcut |
|
| 250 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) & |
|
| 251 |
*(rij/ann(3)-1.0))*fcut*(-2.0*ann(1)/ann(3))*(xj(im)-xi(im))/rij |
|
| 252 |
|
|
| 253 |
df(i,1,im)=df(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*tf & |
|
| 254 |
*(xi(im)-xj(im))/rij |
|
| 255 |
df(i,2,im)=df(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut |
|
| 256 |
df(i,3,im)=df(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*tf & |
|
| 257 |
*(xi(im)-xj(im))/rij |
|
| 258 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) & |
|
| 259 |
*(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 260 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) & |
|
| 261 |
*(rij/ann(3)-1.0))*fcut |
|
| 262 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) & |
|
| 263 |
*(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 264 |
enddo |
|
| 265 |
! |
|
| 266 |
endif |
|
| 267 |
|
|
| 268 |
!**********************************Ni-H interaction |
|
| 269 |
! ********** |
|
| 270 |
if(n1==0.and.n3==0.and.n2==1.and.n4==0)then !Ni-H |
|
| 271 |
call fctNiH(rnh,rij,fcut) |
|
| 272 |
call dfcutNiH(rnh,rij,tf) |
|
| 273 |
eb=eb+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 274 |
er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut |
|
| 275 |
|
|
| 276 |
!Calculation of FORCE |
|
| 277 |
do im=1,3 |
|
| 278 |
dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut & |
|
| 279 |
*(-anh(2)/anh(3))*(xi(im)-xj(im))/rij |
|
| 280 |
dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 281 |
dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut & |
|
| 282 |
*(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij |
|
| 283 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)* & |
|
| 284 |
exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3)) & |
|
| 285 |
*(xj(im)-xi(im))/rij |
|
| 286 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)* & |
|
| 287 |
exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 288 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)* & |
|
| 289 |
exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3)) & |
|
| 290 |
*(xj(im)-xi(im))/rij |
|
| 291 |
|
|
| 292 |
df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf & |
|
| 293 |
*(xi(im)-xj(im))/rij |
|
| 294 |
df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 295 |
df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf & |
|
| 296 |
*(xi(im)-xj(im))/rij |
|
| 297 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) & |
|
| 298 |
*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 299 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1) & |
|
| 300 |
*(rij/anh(3)-1.0))*fcut |
|
| 301 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1) & |
|
| 302 |
*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 303 |
enddo |
|
| 304 |
! |
|
| 305 |
! |
|
| 306 |
endif |
|
| 307 |
|
|
| 308 |
!**********************************Ni-C interaction |
|
| 309 |
! ********** |
|
| 310 |
if(n1==0.and.n3==0.and.n2==1.and.n4==1)then !Ni-C |
|
| 311 |
call fctNiC(rnc,rij,fcut) |
|
| 312 |
call dfcutNiC(rnc,rij,tf) |
|
| 313 |
eb=eb+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 314 |
er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut |
|
| 315 |
|
|
| 316 |
!Calculation of FORCE |
|
| 317 |
do im=1,3 |
|
| 318 |
dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut & |
|
| 319 |
*(-anc(2)/anc(3))*(xi(im)-xj(im))/rij |
|
| 320 |
dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 321 |
dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut & |
|
| 322 |
*(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij |
|
| 323 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* & |
|
| 324 |
anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij |
|
| 325 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* & |
|
| 326 |
anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 327 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* & |
|
| 328 |
anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij |
|
| 329 |
|
|
| 330 |
df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf & |
|
| 331 |
*(xi(im)-xj(im))/rij |
|
| 332 |
df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 333 |
df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf & |
|
| 334 |
*(xi(im)-xj(im))/rij |
|
| 335 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* & |
|
| 336 |
anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 337 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* & |
|
| 338 |
anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 339 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* & |
|
| 340 |
anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 341 |
enddo |
|
| 342 |
! |
|
| 343 |
endif |
|
| 344 |
!!**********************************H-Ni interaction |
|
| 345 |
! ********** |
|
| 346 |
if(n1==1.and.n3==0.and.n2==0.and.n4==0)then ! H-Ni |
|
| 347 |
call fctNiH(rnh,rij,fcut) |
|
| 348 |
call dfcutNiH(rnh,rij,tf) |
|
| 349 |
eb=eb+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 350 |
er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut |
|
| 351 |
|
|
| 352 |
!Calculation of FORCE |
|
| 353 |
do im=1,3 |
|
| 354 |
dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut & |
|
| 355 |
*(-anh(2)/anh(3))*(xi(im)-xj(im))/rij |
|
| 356 |
dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 357 |
dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut & |
|
| 358 |
*(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij |
|
| 359 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) & |
|
| 360 |
*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xj(im)-xi(im))/rij |
|
| 361 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* & |
|
| 362 |
anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 363 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* & |
|
| 364 |
anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xj(im)-xi(im))/rij |
|
| 365 |
|
|
| 366 |
df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf & |
|
| 367 |
*(xi(im)-xj(im))/rij |
|
| 368 |
df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 369 |
df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf & |
|
| 370 |
*(xi(im)-xj(im))/rij |
|
| 371 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0* & |
|
| 372 |
anh(2)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 373 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* & |
|
| 374 |
anh(1)*(rij/anh(3)-1.0))*fcut |
|
| 375 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* & |
|
| 376 |
anh(1)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 377 |
enddo |
|
| 378 |
! |
|
| 379 |
endif |
|
| 380 |
!!**********************************H-H interaction |
|
| 381 |
! ********** |
|
| 382 |
if(n1==1.and.n3==0.and.n2==1.and.n4==0)then !H-H |
|
| 383 |
call fctHH(rhh,rij,fcut) |
|
| 384 |
call dfcutHH(rhh,rij,tf) |
|
| 385 |
eb=eb+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut |
|
| 386 |
er=er+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut |
|
| 387 |
|
|
| 388 |
!Calculation of FORCE |
|
| 389 |
do im=1,3 |
|
| 390 |
dv(i,1,im)=dv(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut & |
|
| 391 |
*(-ahh(2)/ahh(3))*(xi(im)-xj(im))/rij |
|
| 392 |
dv(i,2,im)=dv(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut |
|
| 393 |
dv(i,3,im)=dv(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut & |
|
| 394 |
*(-2.0*ahh(1)/ahh(3))*(xi(im)-xj(im))/rij |
|
| 395 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* & |
|
| 396 |
ahh(2)*(rij/ahh(3)-1.0))*fcut*(-ahh(2)/ahh(3))*(xj(im)-xi(im))/rij |
|
| 397 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* & |
|
| 398 |
ahh(1)*(rij/ahh(3)-1.0))*fcut |
|
| 399 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* & |
|
| 400 |
ahh(1)*(rij/ahh(3)-1.0))*fcut*(-2.0*ahh(1)/ahh(3))*(xj(im)-xi(im))/rij |
|
| 401 |
|
|
| 402 |
df(i,1,im)=df(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*tf & |
|
| 403 |
*(xi(im)-xj(im))/rij |
|
| 404 |
df(i,2,im)=df(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut |
|
| 405 |
df(i,3,im)=df(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*tf & |
|
| 406 |
*(xi(im)-xj(im))/rij |
|
| 407 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* & |
|
| 408 |
ahh(2)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 409 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* & |
|
| 410 |
ahh(1)*(rij/ahh(3)-1.0))*fcut |
|
| 411 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* & |
|
| 412 |
ahh(1)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 413 |
enddo |
|
| 414 |
! |
|
| 415 |
endif |
|
| 416 |
|
|
| 417 |
|
|
| 418 |
!!**********************************H-C interaction |
|
| 419 |
! ********** |
|
| 420 |
if(n1==1.and.n3==0.and.n2==1.and.n4==1)then !H-C |
|
| 421 |
call fctCH(rch,rij,fcut) |
|
| 422 |
call dfcutCH(rch,rij,tf) |
|
| 423 |
eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 424 |
er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut |
|
| 425 |
|
|
| 426 |
!Calculation of FORCE |
|
| 427 |
do im=1,3 |
|
| 428 |
dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut & |
|
| 429 |
*(-ach(2)/ach(3))*(xi(im)-xj(im))/rij |
|
| 430 |
dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 431 |
dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut & |
|
| 432 |
*(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij |
|
| 433 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* & |
|
| 434 |
ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij |
|
| 435 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* & |
|
| 436 |
ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 437 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* & |
|
| 438 |
ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij |
|
| 439 |
|
|
| 440 |
df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf & |
|
| 441 |
*(xi(im)-xj(im))/rij |
|
| 442 |
df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 443 |
df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf & |
|
| 444 |
*(xi(im)-xj(im))/rij |
|
| 445 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* & |
|
| 446 |
ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 447 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* & |
|
| 448 |
ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 449 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* & |
|
| 450 |
ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 451 |
enddo |
|
| 452 |
! |
|
| 453 |
endif |
|
| 454 |
!!**********************************C-Ni interaction |
|
| 455 |
! ********** |
|
| 456 |
if(n1==1.and.n3==1.and.n2==0.and.n4==0)then !C-Ni |
|
| 457 |
call fctNiC(rnc,rij,fcut) |
|
| 458 |
call dfcutNiC(rnc,rij,tf) |
|
| 459 |
eb=eb+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 460 |
er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut |
|
| 461 |
|
|
| 462 |
!Calculation of FORCE |
|
| 463 |
do im=1,3 |
|
| 464 |
dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut & |
|
| 465 |
*(-anc(2)/anc(3))*(xi(im)-xj(im))/rij |
|
| 466 |
dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 467 |
dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut & |
|
| 468 |
*(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij |
|
| 469 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* & |
|
| 470 |
anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij |
|
| 471 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* & |
|
| 472 |
anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 473 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* & |
|
| 474 |
anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij |
|
| 475 |
|
|
| 476 |
df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf & |
|
| 477 |
*(xi(im)-xj(im))/rij |
|
| 478 |
df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 479 |
df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf & |
|
| 480 |
*(xi(im)-xj(im))/rij |
|
| 481 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* & |
|
| 482 |
anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 483 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* & |
|
| 484 |
anc(1)*(rij/anc(3)-1.0))*fcut |
|
| 485 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* & |
|
| 486 |
anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 487 |
enddo |
|
| 488 |
! |
|
| 489 |
|
|
| 490 |
endif |
|
| 491 |
!!**********************************C-H interaction |
|
| 492 |
! ********** |
|
| 493 |
if(n1==1.and.n3==1.and.n2==1.and.n4==0)then !C-H |
|
| 494 |
call fctCH(rch,rij,fcut) |
|
| 495 |
call dfcutCH(rch,rij,tf) |
|
| 496 |
eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 497 |
er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut |
|
| 498 |
|
|
| 499 |
|
|
| 500 |
!Calculation of FORCE |
|
| 501 |
do im=1,3 |
|
| 502 |
dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut & |
|
| 503 |
*(-ach(2)/ach(3))*(xi(im)-xj(im))/rij |
|
| 504 |
dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 505 |
dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut & |
|
| 506 |
*(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij |
|
| 507 |
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* & |
|
| 508 |
ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij |
|
| 509 |
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* & |
|
| 510 |
ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 511 |
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* & |
|
| 512 |
ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij |
|
| 513 |
|
|
| 514 |
df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf & |
|
| 515 |
*(xi(im)-xj(im))/rij |
|
| 516 |
df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 517 |
df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf & |
|
| 518 |
*(xi(im)-xj(im))/rij |
|
| 519 |
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* & |
|
| 520 |
ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 521 |
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* & |
|
| 522 |
ach(1)*(rij/ach(3)-1.0))*fcut |
|
| 523 |
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* & |
|
| 524 |
ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij |
|
| 525 |
enddo |
|
| 526 |
! |
|
| 527 |
|
|
| 528 |
endif |
|
| 529 |
!!***************************ONE end**************** |
|
| 530 |
|
|
| 531 |
enddo |
|
| 532 |
!!******** CALCULATION ENERGY |
|
| 533 |
! |
|
| 534 |
yf=yf+er-sqrt(eb) |
|
| 535 |
!!********* |
|
| 536 |
do im=1,3 |
|
| 537 |
dvdr(i,im)=dvdr(i,im)+dv(i,1,im)-0.5/(dv(i,2,im))**0.5*dv(i,3,im) |
|
| 538 |
dvdr(i,im)=dvdr(i,im)+df(i,1,im)-0.5/(df(i,2,im))**0.5*df(i,3,im) |
|
| 539 |
do j=1,nf(i+mid) |
|
| 540 |
jm=ka(mf(i+mid,j)) |
|
| 541 |
dvdr(jm,im)=dvdr(jm,im)+dv(jm,1,im)-0.5/(dv(i,2,im))**0.5*dv(jm,3,im) |
|
| 542 |
dvdr(jm,im)=dvdr(jm,im)+df(jm,1,im)-0.5/(df(i,2,im))**0.5*df(jm,3,im) |
|
| 543 |
enddo |
|
| 544 |
enddo |
|
| 545 |
! |
|
| 546 |
enddo |
|
| 547 |
!!*****************************TWO ENDS************************** |
|
| 548 |
!***************************************Long-ranged interaction |
|
| 549 |
!do i=nq-4,nq |
|
| 550 |
! zi=pz(i+mid) !-0.3*aa*br(3,3) |
|
| 551 |
! call fct2(zshort,zlong,zi,fcut2) |
|
| 552 |
! call dfcut2(zshort,zlong,zi,tf) |
|
| 553 |
! yf=yf+fcut2*(a(6)-a(7)/zi**2.0) |
|
| 554 |
|
|
| 555 |
! dvdr(i,3)=dvdr(i,3)+fcut2*(a(7)*2.0/zi**3.0)+tf*(a(6)-a(7)/zi**2.0) |
|
| 556 |
|
|
| 557 |
!enddo |
|
| 558 |
!******************************************************************* |
|
| 559 |
v=yf |
|
| 560 |
do im=1,3 |
|
| 561 |
do jm=1,nq |
|
| 562 |
dvdr(jm,im)=dvdr(jm,im) |
|
| 563 |
enddo |
|
| 564 |
enddo |
|
| 565 |
|
|
| 566 |
!PFL: We convert dvdr into Grad |
|
| 567 |
do jm=1,nq |
|
| 568 |
Iat=Jm |
|
| 569 |
if (renum) Iat=order(jm) |
|
| 570 |
Grad(3*Iat-2:3*Iat)=dvdr(jm,1:3) |
|
| 571 |
enddo |
|
| 572 |
|
|
| 573 |
if (debug) WRITE(*,*) "EGRAD_CHAMFRE, E=",V |
|
| 574 |
|
|
| 575 |
if (debug) call Header('Egrad_chamfre OVER')
|
|
| 576 |
|
|
| 577 |
!********************************************ENDENDENDENDENDENDENDEND |
|
| 578 |
return |
|
| 579 |
end subroutine egrad_chamfre |
|
| 580 |
|
|
| 581 |
!****************************************** |
|
| 582 |
!**********PARAMETERS********************** |
|
| 583 |
subroutine fctNiNi(r2,r,fcut) |
|
| 584 |
|
|
| 585 |
IMPLICIT NONE |
|
| 586 |
! integer, parameter :: KINT = kind(1) |
|
| 587 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 588 |
|
|
| 589 |
REAL(KREAL) :: r2,r,fcut,r1,pi |
|
| 590 |
pi=4.0*atan(1.0) |
|
| 591 |
r1=r2-0.5 |
|
| 592 |
|
|
| 593 |
if(r<=r1)then |
|
| 594 |
fcut=1.0 |
|
| 595 |
endif |
|
| 596 |
if(r>r1.and.r<=r2)then |
|
| 597 |
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 598 |
endif |
|
| 599 |
if(r>r2)then |
|
| 600 |
fcut=0.0 |
|
| 601 |
endif |
|
| 602 |
return |
|
| 603 |
end subroutine fctNiNi |
|
| 604 |
!!************ |
|
| 605 |
subroutine dfcutNiNi(r2,r,df) |
|
| 606 |
IMPLICIT NONE |
|
| 607 |
! integer, parameter :: KINT = kind(1) |
|
| 608 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 609 |
REAL(KREAL) :: r, r2,df,pi,r1 |
|
| 610 |
pi=4.0*atan(1.0) |
|
| 611 |
r1=r2-0.5 |
|
| 612 |
|
|
| 613 |
if(r<=r1)then |
|
| 614 |
df=0.0 |
|
| 615 |
endif |
|
| 616 |
if(r>r1.and.r<=r2)then |
|
| 617 |
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 618 |
endif |
|
| 619 |
if(r>r2)then |
|
| 620 |
df=0.0 |
|
| 621 |
endif |
|
| 622 |
return |
|
| 623 |
end subroutine dfcutNiNi |
|
| 624 |
|
|
| 625 |
|
|
| 626 |
!!************ |
|
| 627 |
subroutine fctNiH(r2,r,fcut) |
|
| 628 |
IMPLICIT NONE |
|
| 629 |
! integer, parameter :: KINT = kind(1) |
|
| 630 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 631 |
|
|
| 632 |
REAL(KREAL) ::r2,r,fcut |
|
| 633 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 634 |
|
|
| 635 |
rshort=2.7 |
|
| 636 |
rlong=3.0 |
|
| 637 |
pi=4.0*atan(1.0) |
|
| 638 |
r1=rshort |
|
| 639 |
r2=rlong |
|
| 640 |
if(r<=r1)then |
|
| 641 |
fcut=1.0 |
|
| 642 |
endif |
|
| 643 |
if(r>r1.and.r<=r2)then |
|
| 644 |
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 645 |
endif |
|
| 646 |
if(r>r2)then |
|
| 647 |
fcut=0.0 |
|
| 648 |
endif |
|
| 649 |
return |
|
| 650 |
end subroutine fctNiH |
|
| 651 |
!!********** |
|
| 652 |
|
|
| 653 |
subroutine dfcutNiH(r2,r,df) |
|
| 654 |
|
|
| 655 |
|
|
| 656 |
IMPLICIT NONE |
|
| 657 |
! integer, parameter :: KINT = kind(1) |
|
| 658 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 659 |
|
|
| 660 |
|
|
| 661 |
REAL(KREAL) ::r2,r,df |
|
| 662 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 663 |
|
|
| 664 |
rshort=2.7 |
|
| 665 |
rlong=3.0 |
|
| 666 |
pi=4.0*atan(1.0) |
|
| 667 |
r1=rshort |
|
| 668 |
r2=rlong |
|
| 669 |
if(r<=r1)then |
|
| 670 |
df=0.0 |
|
| 671 |
endif |
|
| 672 |
if(r>r1.and.r<=r2)then |
|
| 673 |
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 674 |
endif |
|
| 675 |
if(r>r2)then |
|
| 676 |
df=0.0 |
|
| 677 |
endif |
|
| 678 |
return |
|
| 679 |
end subroutine dfcutNiH |
|
| 680 |
|
|
| 681 |
!!************ |
|
| 682 |
subroutine fctNiC(r2,r,fcut) |
|
| 683 |
|
|
| 684 |
IMPLICIT NONE |
|
| 685 |
! integer, parameter :: KINT = kind(1) |
|
| 686 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 687 |
|
|
| 688 |
|
|
| 689 |
REAL(KREAL) ::r2,r,fcut |
|
| 690 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 691 |
|
|
| 692 |
rshort=3.2 |
|
| 693 |
rlong=3.5 |
|
| 694 |
pi=4.0*atan(1.0) |
|
| 695 |
r1=rshort |
|
| 696 |
r2=rlong |
|
| 697 |
if(r<=r1)then |
|
| 698 |
fcut=1.0 |
|
| 699 |
endif |
|
| 700 |
if(r>r1.and.r<=r2)then |
|
| 701 |
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 702 |
endif |
|
| 703 |
if(r>r2)then |
|
| 704 |
fcut=0.0 |
|
| 705 |
endif |
|
| 706 |
return |
|
| 707 |
end subroutine fctNiC |
|
| 708 |
|
|
| 709 |
!!*********** |
|
| 710 |
subroutine dfcutNiC(r2,r,df) |
|
| 711 |
|
|
| 712 |
IMPLICIT NONE |
|
| 713 |
! integer, parameter :: KINT = kind(1) |
|
| 714 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 715 |
|
|
| 716 |
|
|
| 717 |
REAL(KREAL) ::r2,r,df |
|
| 718 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 719 |
|
|
| 720 |
rshort=3.2 |
|
| 721 |
rlong=3.5 |
|
| 722 |
pi=4.0*atan(1.0) |
|
| 723 |
r1=rshort |
|
| 724 |
r2=rlong |
|
| 725 |
if(r<=r1)then |
|
| 726 |
df=0.0 |
|
| 727 |
endif |
|
| 728 |
if(r>r1.and.r<=r2)then |
|
| 729 |
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 730 |
endif |
|
| 731 |
if(r>r2)then |
|
| 732 |
df=0.0 |
|
| 733 |
endif |
|
| 734 |
return |
|
| 735 |
end subroutine dfcutNiC |
|
| 736 |
|
|
| 737 |
!!************ |
|
| 738 |
subroutine fctCH(r2,r,fcut) |
|
| 739 |
|
|
| 740 |
|
|
| 741 |
IMPLICIT NONE |
|
| 742 |
! integer, parameter :: KINT = kind(1) |
|
| 743 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 744 |
|
|
| 745 |
|
|
| 746 |
REAL(KREAL) ::r2,r,fcut |
|
| 747 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 748 |
|
|
| 749 |
rshort=2.1 |
|
| 750 |
rlong=3.2 |
|
| 751 |
pi=4.0*atan(1.0) |
|
| 752 |
r1=rshort |
|
| 753 |
r2=rlong |
|
| 754 |
if(r<=r1)then |
|
| 755 |
fcut=1.0 |
|
| 756 |
endif |
|
| 757 |
if(r>r1.and.r<=r2)then |
|
| 758 |
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 759 |
endif |
|
| 760 |
if(r>r2)then |
|
| 761 |
fcut=0.0 |
|
| 762 |
endif |
|
| 763 |
return |
|
| 764 |
end subroutine fctCH |
|
| 765 |
!!***************** |
|
| 766 |
|
|
| 767 |
subroutine dfcutCH(r2,r,df) |
|
| 768 |
|
|
| 769 |
|
|
| 770 |
IMPLICIT NONE |
|
| 771 |
! integer, parameter :: KINT = kind(1) |
|
| 772 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 773 |
|
|
| 774 |
|
|
| 775 |
REAL(KREAL) ::r2,r,df |
|
| 776 |
REAL(KREAL) ::rshort,rlong,r1,pi |
|
| 777 |
|
|
| 778 |
rshort=2.1 |
|
| 779 |
rlong=3.2 |
|
| 780 |
pi=4.0*atan(1.0) |
|
| 781 |
r1=rshort |
|
| 782 |
r2=rlong |
|
| 783 |
if(r<=r1)then |
|
| 784 |
df=0.0 |
|
| 785 |
endif |
|
| 786 |
if(r>r1.and.r<=r2)then |
|
| 787 |
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 788 |
endif |
|
| 789 |
if(r>r2)then |
|
| 790 |
df=0.0 |
|
| 791 |
endif |
|
| 792 |
return |
|
| 793 |
end subroutine dfcutCH |
|
| 794 |
!!************ |
|
| 795 |
subroutine fctHH(r2,r,fcut) |
|
| 796 |
|
|
| 797 |
IMPLICIT NONE |
|
| 798 |
! integer, parameter :: KINT = kind(1) |
|
| 799 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 800 |
|
|
| 801 |
|
|
| 802 |
REAL(KREAL) ::r2,r,fcut,r1,pi |
|
| 803 |
pi=4.0*atan(1.0) |
|
| 804 |
r1=r2-1.1 |
|
| 805 |
if(r<=r1)then |
|
| 806 |
fcut=1.0 |
|
| 807 |
endif |
|
| 808 |
if(r>r1.and.r<=r2)then |
|
| 809 |
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 810 |
endif |
|
| 811 |
if(r>r2)then |
|
| 812 |
fcut=0.0 |
|
| 813 |
endif |
|
| 814 |
return |
|
| 815 |
end subroutine fctHH |
|
| 816 |
!!************ |
|
| 817 |
subroutine dfcutHH(r2,r,df) |
|
| 818 |
|
|
| 819 |
|
|
| 820 |
IMPLICIT NONE |
|
| 821 |
! integer, parameter :: KINT = kind(1) |
|
| 822 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 823 |
|
|
| 824 |
|
|
| 825 |
REAL(KREAL) ::r2,r,df,r1,pi |
|
| 826 |
pi=4.0*atan(1.0) |
|
| 827 |
r1=r2-1.1 |
|
| 828 |
if(r<=r1)then |
|
| 829 |
df=0.0 |
|
| 830 |
endif |
|
| 831 |
if(r>r1.and.r<=r2)then |
|
| 832 |
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 833 |
endif |
|
| 834 |
if(r>r2)then |
|
| 835 |
df=0.0 |
|
| 836 |
endif |
|
| 837 |
return |
|
| 838 |
end subroutine dfcutHH |
|
| 839 |
|
|
| 840 |
!!*********************LONG INTERACTION |
|
| 841 |
|
|
| 842 |
subroutine dfcut2(zshort,zlong,r,df) |
|
| 843 |
|
|
| 844 |
|
|
| 845 |
IMPLICIT NONE |
|
| 846 |
! integer, parameter :: KINT = kind(1) |
|
| 847 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 848 |
|
|
| 849 |
|
|
| 850 |
REAL(KREAL) ::zshort,zlong,r,df,pi, r1, r2 |
|
| 851 |
|
|
| 852 |
pi=4.0*atan(1.0) |
|
| 853 |
|
|
| 854 |
r1=zshort |
|
| 855 |
r2=zlong |
|
| 856 |
|
|
| 857 |
if(r<=r1)then |
|
| 858 |
df=0.0 |
|
| 859 |
endif |
|
| 860 |
if(r>r1.and.r<=r2)then |
|
| 861 |
df=pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1)) |
|
| 862 |
endif |
|
| 863 |
if(r>r2)then |
|
| 864 |
df=0.0 |
|
| 865 |
endif |
|
| 866 |
|
|
| 867 |
return |
|
| 868 |
end subroutine dfcut2 |
|
| 869 |
|
|
| 870 |
subroutine fct2(zshort,zlong,r,fcut) |
|
| 871 |
|
|
| 872 |
IMPLICIT NONE |
|
| 873 |
! integer, parameter :: KINT = kind(1) |
|
| 874 |
integer, parameter :: KREAL = kind(1.0d0) |
|
| 875 |
|
|
| 876 |
|
|
| 877 |
REAL(KREAL) ::zshort,zlong,r,fcut,pi,r1,r2 |
|
| 878 |
|
|
| 879 |
pi=4.0*atan(1.0) |
|
| 880 |
|
|
| 881 |
r1=zshort |
|
| 882 |
r2=zlong |
|
| 883 |
|
|
| 884 |
if(r<=r1)then |
|
| 885 |
fcut=0.0 |
|
| 886 |
endif |
|
| 887 |
if(r>r1.and.r<=r2)then |
|
| 888 |
fcut=1.0-(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0 |
|
| 889 |
endif |
|
| 890 |
if(r>r2)then |
|
| 891 |
fcut=1.0 |
|
| 892 |
endif |
|
| 893 |
|
|
| 894 |
return |
|
| 895 |
end subroutine fct2 |
|
| 896 |
|
|
| src/Tensor.f90 (revision 12) | ||
|---|---|---|
| 1 |
subroutine tensor(natom,rx,ry,rz,Tinert,CoordG,rmass) |
|
| 2 |
! c |
|
| 3 |
! c elements of inertial tensor for a set of atomic systems. |
|
| 4 |
! c |
|
| 5 |
! c rx,ry,rz : arrays natom with cartesian coordinates of natom particles; input |
|
| 6 |
! c txx,tyy,tzz,txy,txz,tyz: elements of inertial tensors ; output |
|
| 7 |
! c a,b,c : returning center of mass coordinates. |
|
| 8 |
! c rmass : vector of length natom with particle masses; input |
|
| 9 |
! c natom : number of particles per system; input |
|
| 10 |
! c IOOUT : unit number for messages, silent for IOOUT.le.0; input |
|
| 11 |
! c subroutines called: cmshft m. lewerenz jul/92 |
|
| 12 |
! c |
|
| 13 |
|
|
| 14 |
use Io_module |
|
| 15 |
|
|
| 16 |
IMPLICIT NONE |
|
| 17 |
|
|
| 18 |
REAL(KREAL), parameter :: zero=0.d0 |
|
| 19 |
|
|
| 20 |
INTEGER(KINT), INTENT(IN) :: Natom |
|
| 21 |
real(KREAL) :: rx(Natom),ry(Natom),rz(Natom),rmass(Natom) |
|
| 22 |
real(KREAL) :: a,b,c,Tinert(3,3),CoordG(3) |
|
| 23 |
|
|
| 24 |
INTEGER(KINT) :: jover, J |
|
| 25 |
REAL(KREAL) :: Txx, Tyy, Tzz, Txy, Txz, Tyz, Txxi |
|
| 26 |
REAL(KREAL) :: Ryb, Rzc, Rxa, Rmj2, ryb2, rzc2, rxa2, Rmj |
|
| 27 |
|
|
| 28 |
|
|
| 29 |
if(natom.le.0) then |
|
| 30 |
if(IOOUT.gt.0) write(IOOUT,'(/2a/)') & |
|
| 31 |
' *** error, illegal array dimension(s) in tensor,', & |
|
| 32 |
' return without action ***' |
|
| 33 |
else if(natom.eq.1) then |
|
| 34 |
a=rx(1) |
|
| 35 |
b=ry(1) |
|
| 36 |
c=rz(1) |
|
| 37 |
txx=zero |
|
| 38 |
tyy=zero |
|
| 39 |
tzz=zero |
|
| 40 |
txy=zero |
|
| 41 |
txz=zero |
|
| 42 |
tyz=zero |
|
| 43 |
else |
|
| 44 |
!c |
|
| 45 |
!c first center of mass coordinates, then tensor elements |
|
| 46 |
!c |
|
| 47 |
call cmshft(natom,rx,ry,rz,rmass,a,b,c,0) |
|
| 48 |
CoordG(1)=a |
|
| 49 |
CoordG(2)=b |
|
| 50 |
CoordG(3)=c |
|
| 51 |
|
|
| 52 |
jover=mod(natom,2) |
|
| 53 |
if(jover.eq.0) then |
|
| 54 |
txx=zero |
|
| 55 |
tyy=zero |
|
| 56 |
tzz=zero |
|
| 57 |
txy=zero |
|
| 58 |
txz=zero |
|
| 59 |
tyz=zero |
|
| 60 |
else |
|
| 61 |
rmj=rmass(1) |
|
| 62 |
ryb=b-ry(1) |
|
| 63 |
rzc=c-rz(1) |
|
| 64 |
tyz=-rmj*ryb*rzc |
|
| 65 |
rxa=a-rx(1) |
|
| 66 |
txz=(-rmj*rxa)*rzc |
|
| 67 |
txy=(-rmj*rxa)*ryb |
|
| 68 |
txx=rmj*((ryb*ryb)+(rzc*rzc)) |
|
| 69 |
tyy=rmj*((rxa*rxa)+(rzc*rzc)) |
|
| 70 |
tzz=rmj*((rxa*rxa)+(ryb*ryb)) |
|
| 71 |
end if |
|
| 72 |
!c |
|
| 73 |
!c force vectorization of inner loop if necessary (ibm-3090) |
|
| 74 |
!c |
|
| 75 |
do j=jover+1,natom,2 |
|
| 76 |
rmj=rmass(j) |
|
| 77 |
rmj2=rmass(j+1) |
|
| 78 |
ryb=b-ry(j) |
|
| 79 |
rzc=c-rz(j) |
|
| 80 |
txxi=txx+rmj*((ryb*ryb)+(rzc*rzc)) |
|
| 81 |
ryb2=b-ry(j+1) |
|
| 82 |
rzc2=c-rz(j+1) |
|
| 83 |
txx=txxi+rmj2*((ryb2*ryb2)+(rzc2*rzc2)) |
|
| 84 |
tyz=tyz-rmj*ryb*rzc-rmj2*ryb2*rzc2 |
|
| 85 |
rxa=a-rx(j) |
|
| 86 |
rxa2=a-rx(j+1) |
|
| 87 |
txy=txy-(rmj*rxa)*ryb-(rmj2*rxa2)*ryb2 |
|
| 88 |
txz=txz-(rmj*rxa)*rzc-(rmj2*rxa2)*rzc2 |
|
| 89 |
tyy=tyy+rmj*((rxa*rxa)+(rzc*rzc)) & |
|
| 90 |
+rmj2*((rxa2*rxa2)+(rzc2*rzc2)) |
|
| 91 |
tzz=tzz+rmj*((rxa*rxa)+(ryb*ryb)) & |
|
| 92 |
+rmj2*((rxa2*rxa2)+(ryb2*ryb2)) |
|
| 93 |
END DO |
|
| 94 |
end if |
|
| 95 |
Tinert(1,1)=txx |
|
| 96 |
Tinert(1,2)=txy |
|
| 97 |
Tinert(1,3)=txz |
|
| 98 |
Tinert(2,1)=txy |
|
| 99 |
Tinert(2,2)=tyy |
|
| 100 |
Tinert(2,3)=tyz |
|
| 101 |
Tinert(3,1)=txz |
|
| 102 |
Tinert(3,2)=tyz |
|
| 103 |
Tinert(3,3)=tzz |
|
| 104 |
return |
|
| 105 |
end subroutine tensor |
|
| src/Ginvse.f90 (revision 12) | ||
|---|---|---|
| 1 |
|
|
| 2 |
SUBROUTINE GINVSE(A, MA, NA, AIN, NAIN, M, N, TOL, V, NV, FAIL, NOUT) |
|
| 3 |
|
|
| 4 |
Use Io_module |
|
| 5 |
! INVERSE.OF.A = V * INVERSE.OF.S * TRANSPOSE.OF.U |
|
| 6 |
! VIA SINGULAR VALUE DECOMPOSITION USING ALGORITHMS 1 AND 2 OF |
|
| 7 |
! NASH J C, COMPACT NUMERICAL METHODS FOR COMPUTERS, 1979, |
|
| 8 |
! ADAM HILGER, BRISTOL OR HALSTED PRESS, N.Y. |
|
| 9 |
! |
|
| 10 |
! |
|
| 11 |
! CALLING SEQUENCE |
|
| 12 |
! |
|
| 13 |
! CALL GINVSE (A,MA,NA,AIN,NAIN,M,N,TOL,V,NV,FAIL,NOUT) |
|
| 14 |
! |
|
| 15 |
! A = MATRIX OF SIZE M BY N TO BE 'INVERTED' |
|
| 16 |
! A IS DESTROYED BY THIS ROUTINE AND BECOMES THE |
|
| 17 |
! U MATRIX OF THE SINGULAR VALUE DECOMPOSITION |
|
| 18 |
! MA = ROW DIMENSION OF ARRAY A |
|
| 19 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 20 |
! NA = COLUMN DIMENSION OF ARRAY A |
|
| 21 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 22 |
! AIN = COMPUTED 'INVERSE' (N ROWS BY M COLS.) |
|
| 23 |
! NAIN = ROW DIMENSION OF ARRAY AIN |
|
| 24 |
! ( COLUMN DIMENSION ASSUMED TO BE AT LEAST M ) |
|
| 25 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 26 |
! M = NUMBER OF ROWS IN MATRIX A AND |
|
| 27 |
! THE NUMBER OF COLUMNS IN MATRIX AIN |
|
| 28 |
! N = NUMBER OF COLUMNS IN MATRIX A AND |
|
| 29 |
! THE NUMBER OF ROWS IN MATRIX AIN |
|
| 30 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 31 |
! TOL = MACHINE PRECISION FOR COMPUTING ENVIRONMENT |
|
| 32 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 33 |
! V = WORK MATRIX WHICH BECOMES THE V MATRIX (N BY N) |
|
| 34 |
! OF THE SINGULAR VALUE DECOMPOSITION |
|
| 35 |
! NV = DIMENSIONS OF V (BOTH ROWS AND COLUMNS) |
|
| 36 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 37 |
! FAIL = ERROR FLAG - .TRUE. IMPLIES FAILURE OF GINVSE |
|
| 38 |
! NOUT = OUTPUT CHANNEL NUMBER |
|
| 39 |
! UNCHANGED BY THIS SUB-PROGRAM |
|
| 40 |
! |
|
| 41 |
! NOTE. THIS ROUTINE IS NOT DESIGNED TO PRODUCE A 'GOOD' |
|
| 42 |
! GENERALIZED INVERSE. IT IS MEANT ONLY TO FURNISH TEST |
|
| 43 |
! DATA FOR THE ROUTINES GMATX, ZIELKE, AND PTST |
|
| 44 |
! |
|
| 45 |
INTEGER(KINT), INTENT(IN) :: MA, NA, NAIN,N, M,NV, NOUT |
|
| 46 |
INTEGER(KINT) :: ICOUNT, NM1, I, J, JP1, K,SWEEP, LIMIT |
|
| 47 |
LOGICAL, INTENT(INOUT) :: FAIL |
|
| 48 |
REAL(KREAL), INTENT(IN) :: TOL |
|
| 49 |
REAL(KREAL) :: P, Q, R, VV, C, S, TEMP |
|
| 50 |
REAL(KREAL), INTENT(INOUT) :: V(NV,NV) |
|
| 51 |
REAL(KREAL), INTENT(INOUT) :: AIN(NAIN,M) |
|
| 52 |
REAL(KREAL), INTENT(INOUT) :: A(MA,NA) |
|
| 53 |
|
|
| 54 |
LOGICAL :: debug |
|
| 55 |
INTERFACE |
|
| 56 |
function valid(string) result (isValid) |
|
| 57 |
CHARACTER(*), intent(in) :: string |
|
| 58 |
logical :: isValid |
|
| 59 |
END function VALID |
|
| 60 |
|
|
| 61 |
END INTERFACE |
|
| 62 |
|
|
| 63 |
debug=valid('Ginvse')
|
|
| 64 |
if (debug) WRITE(*,*) '=======Entering Ginvse=======' |
|
| 65 |
|
|
| 66 |
!Print *, 'A=' |
|
| 67 |
DO J=1,MA |
|
| 68 |
!WRITE(*,'(50(1X,F10.2))') AIN(:,J) |
|
| 69 |
!Print *, A(:,J) |
|
| 70 |
END DO |
|
| 71 |
!Print *, 'End of A' |
|
| 72 |
! |
|
| 73 |
! INITIALLY SET FAIL=.FALSE. TO IMPLY CORRECT EXECUTION |
|
| 74 |
! |
|
| 75 |
FAIL = .FALSE. |
|
| 76 |
! |
|
| 77 |
! TESTS ON VALIDITY OF DIMENSIONS |
|
| 78 |
! |
|
| 79 |
IF (M.LT.2 .OR. MA.LT.2 .OR. N.LT.2 .OR. NA.LT.2) GO TO 230 |
|
| 80 |
IF (N.GT.NA) GO TO 210 |
|
| 81 |
IF (M.GT.MA) GO TO 220 |
|
| 82 |
! |
|
| 83 |
! N MUST NOT EXCEED M, OTHERWISE GINVSE WILL FAIL |
|
| 84 |
! |
|
| 85 |
IF (N.GT.M) GO TO 200 |
|
| 86 |
! |
|
| 87 |
! SWEEP COUNTER INITIALIZED TO ZERO |
|
| 88 |
! |
|
| 89 |
SWEEP = 0 |
|
| 90 |
! |
|
| 91 |
! SET SWEEP LIMIT |
|
| 92 |
! |
|
| 93 |
LIMIT = MAX0(N,6) |
|
| 94 |
! |
|
| 95 |
! MAX0(N,6) WAS CHOSEN FROM EXPERIENCE |
|
| 96 |
! |
|
| 97 |
! |
|
| 98 |
! V(I,J) INITIALLY N BY N IDENTITY |
|
| 99 |
! |
|
| 100 |
DO 20 I=1,N |
|
| 101 |
DO 10 J=1,N |
|
| 102 |
V(I,J) = 0.0 |
|
| 103 |
10 CONTINUE |
|
| 104 |
V(I,I) = 1.0 |
|
| 105 |
20 CONTINUE |
|
| 106 |
! |
|
| 107 |
! INITIALIZE ROTATION COUNTER (COUNTS DOWN TO 0) |
|
| 108 |
! |
|
| 109 |
30 ICOUNT = N*(N-1)/2 |
|
| 110 |
! |
|
| 111 |
! COUNT SWEEP |
|
| 112 |
! |
|
| 113 |
SWEEP = SWEEP + 1 |
|
| 114 |
NM1 = N - 1 |
|
| 115 |
DO 110 J=1,NM1 |
|
| 116 |
JP1 = J + 1 |
|
| 117 |
DO 100 K=JP1,N |
|
| 118 |
P = 0. |
|
| 119 |
Q = 0. |
|
| 120 |
R = 0. |
|
| 121 |
DO 40 I=1,M |
|
| 122 |
! |
|
| 123 |
! TEST FOR AND AVOID UNDERFLOW |
|
| 124 |
! NOT NEEDED FOR MACHINES WHICH UNDERFLOW TO ZERO WITHOUT |
|
| 125 |
! ERROR MESSAGE |
|
| 126 |
! |
|
| 127 |
IF (ABS(A(I,J)).GE.TOL) Q = Q + A(I,J)*A(I,J) |
|
| 128 |
!Print *, '1 Q=',Q |
|
| 129 |
!Print *, 'A(',I,J,')=',A(I,J)
|
|
| 130 |
IF (ABS(A(I,K)).GE.TOL) R = R + A(I,K)*A(I,K) |
|
| 131 |
IF (ABS(A(I,J)/TOL)*ABS(A(I,K)/TOL).GE.1.0) P = P + A(I,J)*A(I,K) |
|
| 132 |
40 CONTINUE |
|
| 133 |
!Print *, '2 Q=',Q |
|
| 134 |
IF (Q.GE.R) GO TO 50 |
|
| 135 |
C = 0. |
|
| 136 |
S = 1. |
|
| 137 |
GO TO 60 |
|
| 138 |
50 IF (ABS(Q).LT.TOL .AND. ABS(R).LT.TOL) GO TO 90 |
|
| 139 |
IF (R.EQ.0.0) GO TO 90 |
|
| 140 |
IF ((P/Q)*(P/R).LT.TOL) GO TO 90 |
|
| 141 |
! |
|
| 142 |
! CALCULATE THE SINE AND COSINE OF THE ANGLE OF ROTATION |
|
| 143 |
! |
|
| 144 |
!Print *, '3 Q=',Q |
|
| 145 |
!Print *, 'R=',R |
|
| 146 |
Q = Q - R |
|
| 147 |
VV = SQRT(4.*P*P+Q*Q) |
|
| 148 |
!Print *, 'VV=',VV |
|
| 149 |
!Print *, 'P=',P |
|
| 150 |
!Print *, 'Q=',Q |
|
| 151 |
C = SQRT((VV+Q)/(2.*VV)) |
|
| 152 |
!Print *, 'C=', C |
|
| 153 |
S = P/(VV*C) |
|
| 154 |
!Print *, 'S=', S |
|
| 155 |
! |
|
| 156 |
! APPLY THE ROTATION TO A |
|
| 157 |
! |
|
| 158 |
60 DO 70 I=1,M |
|
| 159 |
R = A(I,J) |
|
| 160 |
A(I,J) = R*C + A(I,K)*S |
|
| 161 |
A(I,K) = -R*S + A(I,K)*C |
|
| 162 |
70 CONTINUE |
|
| 163 |
!if (debug) WRITE(*,*) '=======After rotation of A=======' |
|
| 164 |
! |
|
| 165 |
! APPLY THE ROTATION TO V |
|
| 166 |
! |
|
| 167 |
DO 80 I=1,N |
|
| 168 |
R = V(I,J) |
|
| 169 |
V(I,J) = R*C + V(I,K)*S |
|
| 170 |
!Print *, 'V(I,J)=', V(I,J) |
|
| 171 |
!Print *, 'R=', R |
|
| 172 |
V(I,K) = -R*S + V(I,K)*C |
|
| 173 |
!Print *, 'V(I,K)=', V(I,K) |
|
| 174 |
80 CONTINUE |
|
| 175 |
GO TO 100 |
|
| 176 |
90 ICOUNT = ICOUNT - 1 |
|
| 177 |
100 CONTINUE |
|
| 178 |
110 CONTINUE |
|
| 179 |
!if (debug) WRITE(*,*) '=======After rotation of V=======' |
|
| 180 |
! |
|
| 181 |
! PRINT THE NUMBER OF SWEEPS AND ROTATIONS |
|
| 182 |
! |
|
| 183 |
IF (NOUT.GT.0) WRITE (NOUT,9997) SWEEP, ICOUNT |
|
| 184 |
! |
|
| 185 |
! CHECK NUMBER OF SWEEPS AND ROTATIONS (TERMINATION TEST) |
|
| 186 |
! |
|
| 187 |
IF (ICOUNT.GT.0 .AND. SWEEP.LT.LIMIT) GO TO 30 |
|
| 188 |
IF (SWEEP.GE.LIMIT .AND. NOUT.GT.0) WRITE (NOUT,9996) |
|
| 189 |
DO 160 J=1,N |
|
| 190 |
Q = 0. |
|
| 191 |
DO 120 I=1,M |
|
| 192 |
Q = Q + A(I,J)*A(I,J) |
|
| 193 |
120 CONTINUE |
|
| 194 |
! |
|
| 195 |
! ARBITRARY RANK DECISION |
|
| 196 |
! |
|
| 197 |
IF (J.EQ.1) VV = 1.0E-3*SQRT(Q) |
|
| 198 |
IF (SQRT(Q).LE.VV) GO TO 140 |
|
| 199 |
DO 130 I=1,M |
|
| 200 |
A(I,J) = A(I,J)/Q |
|
| 201 |
130 CONTINUE |
|
| 202 |
GO TO 160 |
|
| 203 |
140 DO 150 I=1,M |
|
| 204 |
A(I,J) = 0.0 |
|
| 205 |
150 CONTINUE |
|
| 206 |
160 CONTINUE |
|
| 207 |
|
|
| 208 |
!if (debug) WRITE(*,*) '=======After termination test=======' |
|
| 209 |
|
|
| 210 |
DO 190 I=1,N |
|
| 211 |
DO 180 J=1,M |
|
| 212 |
TEMP = 0.0 |
|
| 213 |
DO 170 K=1,N |
|
| 214 |
TEMP = TEMP + V(I,K)*A(J,K) |
|
| 215 |
!Print *, 'V(I,K)*A(J,K)=', V(I,K)*A(J,K) |
|
| 216 |
!Print *, 'V(I,K)=', V(I,K) |
|
| 217 |
!Print *, 'A(J,K)=', A(J,K) |
|
| 218 |
170 CONTINUE |
|
| 219 |
AIN(I,J) = TEMP |
|
| 220 |
!Print *, 'AIN(I,J)=', AIN(I,J) |
|
| 221 |
180 CONTINUE |
|
| 222 |
190 CONTINUE |
|
| 223 |
RETURN |
|
| 224 |
|
|
| 225 |
!Print *, 'AIN=' |
|
| 226 |
!DO J=1,NAIN |
|
| 227 |
!WRITE(*,'(50(1X,F10.2))') AIN(:,J) |
|
| 228 |
!Print *, AIN(:,J) |
|
| 229 |
!END DO |
|
| 230 |
!Print *, 'End of AIN' |
|
| 231 |
|
|
| 232 |
! |
|
| 233 |
! IF AN ERROR OCCURED THE APPROPRIATE MESSAGE IS PRINTED |
|
| 234 |
! AND FAIL IS SET TO .TRUE. |
|
| 235 |
! |
|
| 236 |
200 IF (NOUT.GT.0) WRITE (NOUT,9995) |
|
| 237 |
FAIL = .TRUE. |
|
| 238 |
RETURN |
|
| 239 |
210 IF (NOUT.GT.0) WRITE (NOUT,9999) N, NA |
|
| 240 |
FAIL = .TRUE. |
|
| 241 |
RETURN |
|
| 242 |
220 IF (NOUT.GT.0) WRITE (NOUT,9998) M, MA |
|
| 243 |
FAIL = .TRUE. |
|
| 244 |
RETURN |
|
| 245 |
230 IF (NOUT.GT.0) WRITE (NOUT,9994) M, MA, N, NA |
|
| 246 |
FAIL = .TRUE. |
|
| 247 |
RETURN |
|
| 248 |
! |
|
| 249 |
! * * * * * * * * * * * * * * * * * * * |
|
| 250 |
! |
|
| 251 |
9999 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR N IN GINVSE NA= , I5, 5(/)) |
|
| 252 |
9998 FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR M IN GINVSE MA= , I5, 5(/)) |
|
| 253 |
9997 FORMAT (8H SWEEP =, I4, 2H , I4, 19H JACOBI ROTATIONS P,8HERFORMED) |
|
| 254 |
9996 FORMAT (21H SWEEP LIMIT REACHED ) |
|
| 255 |
9995 FORMAT (38H0**ERROR** N GREATER THAN M IS INVALID,21H AND GINVSE WILL FAIL, 5(/)) |
|
| 256 |
9994 FORMAT (46H0**ERROR** SOME DIMENSIONS ARE LESS THAN TWO I,8HN GINVSE & |
|
| 257 |
/10X, 3HM =, I5, 5H MA =,I5, 4H N =, I5,5H NA =, I5, 5(/)) |
|
| 258 |
! |
|
| 259 |
! * * * END OF SUBROUTINE GINVSE * * * |
|
| 260 |
! |
|
| 261 |
END SUBROUTINE GINVSE |
|
| 262 |
|
|
| src/refsor.f90 (revision 12) | ||
|---|---|---|
| 1 |
!!! Downloaded from http://www.fortran-2000.com/rank/index.html |
|
| 2 |
!!! on October 13th, 2010 |
|
| 3 |
|
|
| 4 |
!PFL 2010 Oct 13-> |
|
| 5 |
! I deleted the 'module' part to keep only |
|
| 6 |
! the Double subroutine. |
|
| 7 |
|
|
| 8 |
! Module m_refsor |
|
| 9 |
! ! PFL 2010 Oct 13 -> |
|
| 10 |
! !Integer, Parameter :: kdp = selected_real_kind(15) |
|
| 11 |
! Integer, Parameter :: kdp = KIND(1.0D0) |
|
| 12 |
! ! <- PFL 2010 Oct 13 |
|
| 13 |
! public :: refsor |
|
| 14 |
! private :: kdp |
|
| 15 |
! private :: R_refsor, I_refsor, D_refsor |
|
| 16 |
! private :: R_inssor, I_inssor, D_inssor |
|
| 17 |
! private :: R_subsor, I_subsor, D_subsor |
|
| 18 |
! interface refsor |
|
| 19 |
! module procedure d_refsor, r_refsor, i_refsor |
|
| 20 |
! end interface refsor |
|
| 21 |
! contains |
|
| 22 |
|
|
| 23 |
! Subroutine D_refsor (XDONT) |
|
| 24 |
Subroutine Drefsor (XDONT) |
|
| 25 |
! Sorts XDONT into ascending order - Quicksort |
|
| 26 |
! __________________________________________________________ |
|
| 27 |
! Quicksort chooses a "pivot" in the set, and explores the |
|
| 28 |
! array from both ends, looking for a value > pivot with the |
|
| 29 |
! increasing index, for a value <= pivot with the decreasing |
|
| 30 |
! index, and swapping them when it has found one of each. |
|
| 31 |
! The array is then subdivided in 2 ([3]) subsets: |
|
| 32 |
! { values <= pivot} {pivot} {values > pivot}
|
|
| 33 |
! One then call recursively the program to sort each subset. |
|
| 34 |
! When the size of the subarray is small enough, one uses an |
|
| 35 |
! insertion sort that is faster for very small sets. |
|
| 36 |
! Michel Olagnon - Apr. 2000 |
|
| 37 |
! __________________________________________________________ |
|
| 38 |
! __________________________________________________________ |
|
| 39 |
|
|
| 40 |
use VarTypes |
|
| 41 |
|
|
| 42 |
interface |
|
| 43 |
|
|
| 44 |
Recursive Subroutine D_subsor (XDONT, IDEB1, IFIN1) |
|
| 45 |
|
|
| 46 |
use VarTypes |
|
| 47 |
|
|
| 48 |
IMPLICIT NONE |
|
| 49 |
|
|
| 50 |
Real(KREAL), dimension (:), Intent (InOut) :: XDONT |
|
| 51 |
Integer(KINT), Intent (In) :: IDEB1, IFIN1 |
|
| 52 |
|
|
| 53 |
Integer(KINT), Parameter :: NINS = 16 ! Max for insertion sort |
|
| 54 |
Integer(KINT) :: ICRS, IDEB, IDCR, IFIN, IMIL |
|
| 55 |
Real(KREAL) :: XPIV, XWRK |
|
| 56 |
end Subroutine D_subsor |
|
| 57 |
|
|
| 58 |
Subroutine D_inssor (XDONT) |
|
| 59 |
|
|
| 60 |
use VarTypes |
|
| 61 |
|
|
| 62 |
IMPLICIT NONE |
|
| 63 |
Real(KREAL), dimension (:), Intent (InOut) :: XDONT |
|
| 64 |
|
|
| 65 |
Integer(KINT) :: ICRS, IDCR |
|
| 66 |
Real(KREAL) :: XWRK |
|
| 67 |
|
|
| 68 |
End Subroutine D_inssor |
|
| 69 |
|
|
| 70 |
end interface |
|
| 71 |
|
|
| 72 |
Real (KREAL), Dimension (:), Intent (InOut) :: XDONT |
|
| 73 |
! __________________________________________________________ |
|
| 74 |
! |
|
| 75 |
! |
|
| 76 |
Call D_subsor (XDONT, 1, Size (XDONT)) |
|
| 77 |
Call D_inssor (XDONT) |
|
| 78 |
Return |
|
| 79 |
End Subroutine Drefsor |
|
| 80 |
|
|
| 81 |
Recursive Subroutine D_subsor (XDONT, IDEB1, IFIN1) |
|
| 82 |
! Sorts XDONT from IDEB1 to IFIN1 |
|
| 83 |
! __________________________________________________________ |
|
| 84 |
|
|
| 85 |
use VarTypes |
|
| 86 |
|
|
Formats disponibles : Unified diff