subroutine init_dfa_vars include 'DIMENSIONS' include 'COMMON.DFA' integer ii C Number of restraints idisnum = 0 iphinum = 0 ithenum = 0 ineinum = 0 idislis = 0 iphilis = 0 ithelis = 0 ineilis = 0 jneilis = 0 jneinum = 0 kshell = 0 fnei = 0 C For beta nca = 0 icaidx = 0 C real variables CC WEIGHTS for each min sccdist = 0.0d0 fdist = 0.0d0 sccphi = 0.0d0 sccthe = 0.0d0 sccnei = 0.0d0 fphi1 = 0.0d0 fphi2 = 0.0d0 fthe1 = 0.0d0 fthe2 = 0.0d0 C energies edfatot = 0.0d0 edfadis = 0.0d0 edfaphi = 0.0d0 edfathe = 0.0d0 edfanei = 0.0d0 edfabet = 0.0d0 C gradients gdfad = 0.0d0 gdfat = 0.0d0 gdfan = 0.0d0 gdfab = 0.0d0 C weights for each E term C these should be identical with dis_inc = 0.0d0 phi_inc = 0.0d0 the_inc = 0.0d0 nei_inc = 0.0d0 beta_inc = 0.0d0 wshet = 0.0d0 C precalculate exp table! dfaexp = 0.0d0 do ii = 1, 15001 dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0) end do return end subroutine read_dfa_info C C read fragment informations C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.CHAIN' include 'COMMON.DFA' C NOTE THAT FILENAMES are FIXED, CURRENTLY!! C THIS SHOULD BE MODIFIED!! character*320 buffer integer iodfa parameter(iodfa=89) integer i, j, nval integer ica1, ica2,ica3,ica4,ica5 integer ishell, inca, itmp,iitmp double precision wtmp C C READ DISTANCE C open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33) goto 34 33 write(iout,'(a)') 'Error opening dist_dfa.dat file' stop 34 continue write(iout,'(a)') 'dist_dfa.dat is opened!' C read title read(iodfa, '(a)') buffer C read number of restraints read(iodfa, '(i)') IDFADIS read(iodfa, *) dis_inc do i=1, idfadis read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval idisnum(i)=nval idislis(1,i)=ica1 idislis(2,i)=ica2 do j=1, nval read(iodfa,*) tmp fdist(i,j) = tmp enddo do j=1, nval read(iodfa,*) tmp sccdist(i,j) = tmp enddo enddo close(iodfa) write(iout,'(a)') 'dist_dfa.dat is closed!' C READ ANGLE RESTRAINTS C PHI RESTRAINTS open(iodfa, file='phi_dfa.dat',status='old',err=35) goto 36 35 write(iout,'(a)') 'Error opening dist_dfa.dat file' stop 36 continue write(iout,'(a)') 'phi_dfa.dat is opened!' C READ TITLE read(iodfa, '(a)') buffer C READ NUMBER OF RESTRAINTS READ(iodfa, '(i)') IDFAPHI read(iodfa,*) phi_inc do i=1, idfaphi read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval iphinum(i)=nval iphilis(1,i)=ica1 iphilis(2,i)=ica2 iphilis(3,i)=ica3 iphilis(4,i)=ica4 iphilis(5,i)=ica5 do j=1, nval read(iodfa,*) tmp1,tmp2 fphi1(i,j) = tmp1 fphi2(i,j) = tmp2 enddo do j=1, nval read(iodfa,*) tmp sccphi(i,j) = tmp enddo enddo close(iodfa) C THETA RESTRAINTS open(iodfa, file='theta_dfa.dat',status='old',err=41) goto 42 41 write(iout,'(a)') 'Error opening dist_dfa.dat file' stop 42 continue write(iout,'(a)') 'theta_dfa.dat is opened!' C READ TITLE read(iodfa, '(a)') buffer C READ NUMBER OF RESTRAINTS READ(iodfa, '(i)') IDFATHE read(iodfa,*) the_inc do i=1, idfathe read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval ithenum(i)=nval ithelis(1,i)=ica1 ithelis(2,i)=ica2 ithelis(3,i)=ica3 ithelis(4,i)=ica4 ithelis(5,i)=ica5 do j=1, nval read(iodfa,*) tmp1,tmp2 fthe1(i,j) = tmp1 fthe2(i,j) = tmp2 enddo do j=1, nval read(iodfa,*) tmp sccthe(i,j) = tmp enddo enddo close(iodfa) C END of READING ANGLE RESTRAINT! C NUMBER OF NEIGHBOR CAs open(iodfa,file='nei_dfa.dat',status='old',err=37) goto 38 37 write(iout,'(a)') 'Error opening nei_dfa.dat file' stop 38 continue write(iout,'(a)') 'nei_dfa.dat is opened!' C READ TITLE read(iodfa, '(a)') buffer C READ NUMBER OF RESTRAINTS READ(iodfa, '(i)') idfanei read(iodfa,*) nei_inc do i=1, idfanei read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval ineilis(i)=ica1 kshell(i)=ishell ineinum(i)=nval do j=1, nval read(iodfa,*) inca fnei(i,j) = inca C write(*,*) 'READ NEI:',i,j,fnei(i,j) enddo do j=1, nval read(iodfa,*) tmp sccnei(i,j) = tmp enddo enddo close(iodfa) C END OF NEIGHBORING CA C READ BETA RESTRAINT open(iodfa, file='beta_dfa.dat',status='old',err=39) goto 40 39 write(iout,'(a)') 'Error opening beta_dfa.dat file' stop 40 continue write(iout,'(a)') 'beta_dfa.dat is opened!' read(iodfa,'(a)') buffer read(iodfa,'(i)') itmp read(iodfa,*) beta_inc do i=1,itmp read(iodfa,*) ica1, iitmp do j=1,itmp read(iodfa,*) wtmp wshet(i,j) = wtmp c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j) enddo enddo close(iodfa) C END OF BETA RESTRAINT return END subroutine edfad(edfadis) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.DFA' double precision edfadis integer i, iatm1, iatm2, idiff double precision ckk, sckk,dist,texp double precision jix,jiy,jiz,ep,fp,scc gdfad=0.0d0 do i=1, idfadis iatm1=idislis(1,i)+1 iatm2=idislis(2,i)+1 idiff=abs(iatm1-iatm2) jix=c(1,iatm2)-c(1,iatm1) jiy=c(2,iatm2)-c(2,iatm1) jiz=c(3,iatm2)-c(3,iatm1) dist=sqrt(jix*jix+jiy*jiy+jiz*jiz) ckk=ck(idiff) sckk=sck(idiff) scc = 0.0d0 ep = 0.0d0 fp = 0.0d0 write(*,*) 'DIST:', i, dist write(*,*) 'JIX,JIY,JIZ:', jix, jiy, jiz do j = 1, idisnum(i) dd = dist-fdist(i,j) dtmp = dd*dd/ckk if (dtmp.ge.15.0d0) then texp = 0.0d0 else texp = dfaexp( idint(dtmp*1000)+1 )/sckk endif ep=ep+sccdist(i,j)*texp fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk scc=scc+sccdist(i,j) write(*,*) 'i,j:', i,j write(*,*) 'iatm1,iatm2:', iatm1, iatm2 write(*,*) 'dd,fdist:', dd, fdist(i,j) write(*,*) 'ck,sck:', ckk, sckk write(*,*) 'ep,fp:', ep, fp write(*,*) 'scc:', scc enddo ep = -ep/scc*wwdist fp = fp/scc*wwdist if(abs(ep).lt.1.0d-20)then ep=0.0d0 fp=0.0d0 endif write(*,*) '' write(*,*) 'ep,fp:', ep, fp write(*,*) 'edfadis:', edfadis write(*,*) 'gdfad:', gdfad(:,iatm1) write(*,*) 'gdfad:', gdfad(:,iatm2) edfadis = edfadis + ep*dis_inc gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc write(*,*) 'ep,fp:', ep, fp write(*,*) 'edfadis:', edfadis write(*,*) 'gdfad:', gdfad(:,iatm1) write(*,*) 'gdfad:', gdfad(:,iatm2) C gdfad(1,iatm1) = gdfad(1,iatm1)+jix/dist*fp*dis_inc*wwdist C gdfad(2,iatm1) = gdfad(2,iatm1)+jiy/dist*fp*dis_inc*wwdist C gdfad(3,iatm1) = gdfad(3,iatm1)+jiz/dist*fp*dis_inc*wwdist C gdfad(1,iatm2) = gdfad(1,iatm2)-jix/dist*fp*dis_inc*wwdist C gdfad(2,iatm2) = gdfad(2,iatm2)-jiy/dist*fp*dis_inc*wwdist C gdfad(3,iatm2) = gdfad(3,iatm2)-jiz/dist*fp*dis_inc*wwdist enddo print*, 'EDFAD:', edfadis print*, 'COOR' do i=1,nres write(*,'(a,i6,1x,3f12.5)')'DFADC:',i,c(1:3,i) enddo print*, 'FORCE' do i=1,nres write(*,'(a,i6,1x,3f14.7)')'DFADG:', i, gdfad(1:3,i) enddo return end subroutine edfat(edfator) C DFA torsion angle implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.DFA' integer i,j,ii, iii integer iatom(5) double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5) double precision cwidth, cwidth2 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0) edfator= 0.0d0 enephi = 0.0d0 enethe = 0.0d0 gdfat(:,:) = 0.0d0 C START OF PHI ANGLE do i=1, idfaphi aphi = 0.0d0 do iii=1, 5 iatom(iii) = iphilis(iii,i)+1 end do C ANGLE VECTOR CALCULTION RIX=C(1,IATOM(2))-C(1,IATOM(1)) RIY=C(2,IATOM(2))-C(2,IATOM(1)) RIZ=C(3,IATOM(2))-C(3,IATOM(1)) RIPX=C(1,IATOM(3))-C(1,IATOM(2)) RIPY=C(2,IATOM(3))-C(2,IATOM(2)) RIPZ=C(3,IATOM(3))-C(3,IATOM(2)) RIPPX=C(1,IATOM(4))-C(1,IATOM(3)) RIPPY=C(2,IATOM(4))-C(2,IATOM(3)) RIPPZ=C(3,IATOM(4))-C(3,IATOM(3)) RIP3X=C(1,IATOM(5))-C(1,IATOM(4)) RIP3Y=C(2,IATOM(5))-C(2,IATOM(4)) RIP3Z=C(3,IATOM(5))-C(3,IATOM(4)) GIX=RIY*RIPZ-RIZ*RIPY GIY=RIZ*RIPX-RIX*RIPZ GIZ=RIX*RIPY-RIY*RIPX GIPX=RIPY*RIPPZ-RIPZ*RIPPY GIPY=RIPZ*RIPPX-RIPX*RIPPZ GIPZ=RIPX*RIPPY-RIPY*RIPPX CIPX=C(1,IATOM(3))-C(1,IATOM(1)) CIPY=C(2,IATOM(3))-C(2,IATOM(1)) CIPZ=C(3,IATOM(3))-C(3,IATOM(1)) CIPPX=C(1,IATOM(4))-C(1,IATOM(2)) CIPPY=C(2,IATOM(4))-C(2,IATOM(2)) CIPPZ=C(3,IATOM(4))-C(3,IATOM(2)) CIP3X=C(1,IATOM(5))-C(1,IATOM(3)) CIP3Y=C(2,IATOM(5))-C(2,IATOM(3)) CIP3Z=C(3,IATOM(5))-C(3,IATOM(3)) DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ) DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ) DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ) DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z) C END OF ANGLE VECTOR CALCULTION TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ APHI(1)=TDOT/(DGI*DRIPP) TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z APHI(2)=TDOT/(DGIP*DRIP3) ephi = 0.0d0 tfphi1=0.0d0 tfphi2=0.0d0 scc=0.0d0 do j=1, iphinum(i) DDPS1=APHI(1)-FPHI1(i,j) DDPS2=APHI(2)-FPHI2(i,j) DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2 if (dtmp.ge.15.0d0) then ps_tmp = 0.0d0 else ps_tmp = dfaexp(idint(dtmp*1000)+1) endif ephi=ephi+sccphi(i,j)*ps_tmp tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp scc=scc+sccphi(i,j) C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j), C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j) ENDDO ephi=-ephi/scc*phi_inc*wwangle tfphi1=tfphi1/scc*phi_inc*wwangle tfphi2=tfphi2/scc*phi_inc*wwangle IF (ABS(EPHI).LT.1d-20) THEN EPHI=0.0D0 ENDIF IF (ABS(TFPHI1).LT.1d-20) THEN TFPHI1=0.0D0 ENDIF IF (ABS(TFPHI2).LT.1d-20) THEN TFPHI2=0.0D0 ENDIF C FORCE DIRECTION CALCULATION TDX(1:5)=0.0D0 TDY(1:5)=0.0D0 TDZ(1:5)=0.0D0 DM1=1.0d0/(DGI*DRIPP) GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ DM2=GIRPP/(DGI**3*DRIPP) DM3=GIRPP/(DGI*DRIPP**3) DM4=1.0d0/(DGIP*DRIP3) GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z DM5=GIRP3/(DGIP**3*DRIP3) DM6=GIRP3/(DGIP*DRIP3**3) C FIRST ATOM BY PHI1 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1 & +( GIZ* RIPY- GIY* RIPZ)*DM2 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1 & +( GIX* RIPZ- GIZ* RIPX)*DM2 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1 & +( GIY* RIPX- GIX* RIPY)*DM2 TDX(1)=TDX(1)*TFPHI1 TDY(1)=TDY(1)*TFPHI1 TDZ(1)=TDZ(1)*TFPHI1 C SECOND ATOM BY PHI1 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1 & -(CIPY*GIZ-CIPZ*GIY)*DM2 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1 & -(CIPZ*GIX-CIPX*GIZ)*DM2 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1 & -(CIPX*GIY-CIPY*GIX)*DM2 TDX(2)=TDX(2)*TFPHI1 TDY(2)=TDY(2)*TFPHI1 TDZ(2)=TDZ(2)*TFPHI1 C SECOND ATOM BY PHI2 TDX(2)=TDX(2)+ & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2 TDY(2)=TDY(2)+ & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2 TDZ(2)=TDZ(2)+ & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2 C THIRD ATOM BY PHI1 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3 TDX(3)=TDX(3)*TFPHI1 TDY(3)=TDY(3)*TFPHI1 TDZ(3)=TDZ(3)*TFPHI1 C THIRD ATOM BY PHI2 TDX(3)=TDX(3)+ & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2 TDY(3)=TDY(3)+ & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2 TDZ(3)=TDZ(3)+ & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2 C FOURTH ATOM BY PHI1 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1 C FOURTH ATOM BY PHI2 TDX(4)=TDX(4)+ & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5 & + RIP3X*DM6)*TFPHI2 TDY(4)=TDY(4)+ & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5 & + RIP3Y*DM6)*TFPHI2 TDZ(4)=TDZ(4)+ & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4 & -( GIPX*RIPY-RIPX*GIPY)*DM5 & + RIP3Z*DM6)*TFPHI2 C FIFTH ATOM BY PHI2 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2 C END OF FORCE DIRECTION c force calcuation do ii=1,5 gdfat(1,iatom(ii))=gdfat(1,iatom(ii))+tdx(ii) gdfat(2,iatom(ii))=gdfat(2,iatom(ii))+tdy(ii) gdfat(3,iatom(ii))=gdfat(3,iatom(ii))+tdz(ii) enddo c energy calculation enephi = enephi + ephi c end of single assignment statement ENDDO C END OF PHI RESTRAINT C START OF THETA ANGLE do i=1, idfathe athe = 0.0d0 do iii=1, 5 iatom(iii) = iphilis(iii,i)+1 end do c iatom(1:5)=ithelis(1:5,i) C ANGLE VECTOR CALCULTION RIX=C(1,IATOM(2))-C(1,IATOM(1)) RIY=C(2,IATOM(2))-C(2,IATOM(1)) RIZ=C(3,IATOM(2))-C(3,IATOM(1)) RIPX=C(1,IATOM(3))-C(1,IATOM(2)) RIPY=C(2,IATOM(3))-C(2,IATOM(2)) RIPZ=C(3,IATOM(3))-C(3,IATOM(2)) RIPPX=C(1,IATOM(4))-C(1,IATOM(3)) RIPPY=C(2,IATOM(4))-C(2,IATOM(3)) RIPPZ=C(3,IATOM(4))-C(3,IATOM(3)) RIP3X=C(1,IATOM(5))-C(1,IATOM(4)) RIP3Y=C(2,IATOM(5))-C(2,IATOM(4)) RIP3Z=C(3,IATOM(5))-C(3,IATOM(4)) GIX=RIY*RIPZ-RIZ*RIPY GIY=RIZ*RIPX-RIX*RIPZ GIZ=RIX*RIPY-RIY*RIPX GIPX=RIPY*RIPPZ-RIPZ*RIPPY GIPY=RIPZ*RIPPX-RIPX*RIPPZ GIPZ=RIPX*RIPPY-RIPY*RIPPX GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X CIPX=C(1,IATOM(3))-C(1,IATOM(1)) CIPY=C(2,IATOM(3))-C(2,IATOM(1)) CIPZ=C(3,IATOM(3))-C(3,IATOM(1)) CIPPX=C(1,IATOM(4))-C(1,IATOM(2)) CIPPY=C(2,IATOM(4))-C(2,IATOM(2)) CIPPZ=C(3,IATOM(4))-C(3,IATOM(2)) CIP3X=C(1,IATOM(5))-C(1,IATOM(3)) CIP3Y=C(2,IATOM(5))-C(2,IATOM(3)) CIP3Z=C(3,IATOM(5))-C(3,IATOM(3)) DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ) DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ) DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ) DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ) DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z) C END OF ANGLE VECTOR CALCULTION TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ ATHE(1)=TDOT/(DGI*DGIP) TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ ATHE(2)=TDOT/(DGIP*DGIPP) ETHE=0.0D0 TFTHE1=0.0D0 TFTHE2=0.0D0 SCC=0.0D0 TH_TMP=0.0d0 do j=1,ithenum(i) ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref) ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref) dtmp= (ddth1**2+ddth2**2)/cwidth2 if ( dtmp .ge. 15.0d0) then th_tmp = 0.0d0 else th_tmp = dfaexp ( idint(dtmp*1000)+1 ) end if ethe=ethe+sccthe(i,j)*th_tmp tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1) tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2) scc=scc+sccthe(i,j) C write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j), C & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j) enddo ethe=-ethe/scc*the_inc*wwangle tfthe1=tfthe1/scc*the_inc*wwangle tfthe2=tfthe2/scc*the_inc*wwangle if (abs(ethe).lt.tenm20) then ethe=0.0d0 endif if (abs(tfthe1).lt.tenm20) then tfthe1=0.0d0 endif if (abs(tfthe2).lt.tenm20) then tfthe2=0.0d0 endif TDX(1:5)=0.0D0 TDY(1:5)=0.0D0 TDZ(1:5)=0.0D0 DM1=1.0d0/(DGI*DGIP) DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP) DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3) DM4=1.0d0/(DGIP*DGIPP) DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP) DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3) C FIRST ATOM BY THETA1 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1 C SECOND ATOM BY THETA1 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1 & -(CIPY*GIZ-CIPZ*GIY)*DM2 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1 & -(CIPZ*GIX-CIPX*GIZ)*DM2 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1 & -(CIPX*GIY-CIPY*GIX)*DM2 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1 C SECOND ATOM BY THETA2 TDX(2)=TDX(2)+ & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2 TDY(2)=TDY(2)+ & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2 TDZ(2)=TDZ(2)+ & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2 C THIRD ATOM BY THETA1 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1 & -(GIY*RIZ-GIZ*RIY)*DM2 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1 & -(GIZ*RIX-GIX*RIZ)*DM2 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1 & -(GIX*RIY-GIY*RIX)*DM2 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1 C THIRD ATOM BY THETA2 TDX(3)=TDX(3)+ & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2 TDY(3)=TDY(3)+ & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2 TDZ(3)=TDZ(3)+ & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2 C FOURTH ATOM BY THETA1 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1 C FOURTH ATOM BY THETA2 TDX(4)=TDX(4)+ & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2 TDY(4)=TDY(4)+ & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2 TDZ(4)=TDZ(4)+ & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4 & -(GIPX*RIPY-GIPY*RIPX)*DM5 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2 C FIFTH ATOM BY THETA2 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2 C !! END OF FORCE DIRECTION!!!! do ii=1,5 gdfat(1,iatom(ii))=gdfat(1,iatom(ii))+tdx(ii) gdfat(2,iatom(ii))=gdfat(2,iatom(ii))+tdy(ii) gdfat(3,iatom(ii))=gdfat(3,iatom(ii))+tdz(ii) enddo C energy calculation enethe = enethe + ethe ENDDO edfator = enephi + enethe print*, 'EDFAT:',edfator C print*, 'COOR' C do i=1,nres-1 C write(*,'(a,i6,1x,3f12.5)')'DFATC:',i,c(1:3,i) C enddo print*, 'TOR_FORCE' do i=2,nres-1 write(*,'(a,i6,1x,3f14.7)')'DFATG:',i,gdfat(1:3,i) enddo RETURN END subroutine edfan(edfanei) C DFA neighboring CA restraint implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.DFA' integer i,j,imin integer kshnum, n1atom double precision enenei,tmp_n double precision pai,hpai double precision jix,jiy,jiz,ndiff,snorm_nei double precision t2dx(maxres),t2dy(maxres),t2dz(maxres) double precision dr,dr2,half,ntmp parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0) parameter(pai=3.14159265358979323846D0) parameter(hpai=1.5707963267948966D0) parameter(snorm_nei=0.886226925452758D0) edfanei = 0.0d0 enenei = 0.0d0 gdfan = 0.0d0 c print*, 's1:', s1(:) c print*, 's2:', s2(:) do i=1, idfanei kshnum= kshell(i) n1atom= ineilis(i) + 1 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom tmp_n=0.0d0 ftmp=0.0d0 dnei=0.0d0 dist=0.0d0 t1dx=0.0d0 t1dy=0.0d0 t1dz=0.0d0 t2dx=0.0d0 t2dy=0.0d0 t2dz=0.0d0 c do j = 1, nres do j = 2, nres-1 if (n1atom.eq.j) cycle jix=c(1,j)-c(1,n1atom) jiy=c(2,j)-c(2,n1atom) jiz=c(3,j)-c(3,n1atom) dist=sqrt(jix*jix+jiy*jiy+jiz*jiz) C write(*,*) 'N1ATOM,J,DIST:', n1atom, j, dist if(kshnum.ne.1)then if (dist.lt.s1(kshnum).and. & dist.gt.s2(kshnum-1)) then tmp_n=tmp_n+1.0d0 c write(*,*) 'case1:',tmp_n t1dx=t1dx+0.0d0 t1dy=t1dy+0.0d0 t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 elseif(dist.ge.s1(kshnum).and. & dist.le.s2(kshnum)) then dnei=(dist-s1(kshnum))/dr2*pai tmp_n=tmp_n + half*(1+cos(dnei)) c write(*,*) 'case2:',tmp_n ftmp=-pai*sin(dnei)/dr2/dist/2.0d0 c center atom t1dx=t1dx+jix*ftmp t1dy=t1dy+jiy*ftmp t1dz=t1dz+jiz*ftmp c neighbor atoms t2dx(j)=-jix*ftmp t2dy(j)=-jiy*ftmp t2dz(j)=-jiz*ftmp c elseif(dist.ge.s1(kshnum-1).and. & dist.le.s2(kshnum-1)) then dnei=(dist-s1(kshnum-1))/dr2*pai tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei)) c write(*,*) 'case3:',tmp_n ftmp = hpai*sin(dnei)/dr2/dist c center atom t1dx=t1dx+jix*ftmp t1dy=t1dy+jiy*ftmp t1dz=t1dz+jiz*ftmp c neighbor atoms t2dx(j)=-jix*ftmp t2dy(j)=-jiy*ftmp t2dz(j)=-jiz*ftmp endif elseif(kshnum.eq.1) then if(dist.lt.s1(kshnum))then tmp_n=tmp_n+1.0d0 c write(*,*) 'case4:',tmp_n t1dx=t1dx+0.0d0 t1dy=t1dy+0.0d0 t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 elseif(dist.ge.s1(kshnum).and. & dist.le.s2(kshnum))then dnei=(dist-s1(kshnum))/dr2*pai tmp_n=tmp_n + half*(1+cos(dnei)) c write(*,*) 'case5:',tmp_n ftmp = -hpai*sin(dnei)/dr2/dist c center atom t1dx=t1dx+jix*ftmp t1dy=t1dy+jiy*ftmp t1dz=t1dz+jiz*ftmp c neighbor atoms t2dx(j)=-jix*ftmp t2dy(j)=-jiy*ftmp t2dz(j)=-jiz*ftmp endif endif enddo scc=0.0d0 enei=0.0d0 tmp_fnei=0.0d0 ndiff=0.0d0 do imin=1,ineinum(i) ndiff = tmp_n-fnei(i,imin) dtmp = ndiff*ndiff if (dtmp.ge.15.0d0) then ntmp = 0.0d0 else ntmp = dfaexp( idint(dtmp*1000) + 1 ) end if enei = enei + sccnei(i,imin)*ntmp tmp_fnei = tmp_fnei- & sccnei(i,imin)*ntmp*ndiff*2.0d0 scc=scc+sccnei(i,imin) c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n, c & fnei(i,imin),sccnei(i,imin),enei,scc enddo enei=-enei/scc*snorm_nei*nei_inc*wwnei tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei if (abs(enei).lt.1.0d-20)then enei=0.0d0 endif if (abs(tmp_fnei).lt.1.0d-20) then tmp_fnei=0.0d0 endif c force calculation t1dx=t1dx*tmp_fnei t1dy=t1dy*tmp_fnei t1dz=t1dz*tmp_fnei do j = 2, nres-1 t2dx(j)=t2dx(j)*tmp_fnei t2dy(j)=t2dy(j)*tmp_fnei t2dz(j)=t2dz(j)*tmp_fnei enddo gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz do j = 2 , nres-1 gdfan(1,j)=gdfan(1,j)+t2dx(j) gdfan(2,j)=gdfan(2,j)+t2dy(j) gdfan(3,j)=gdfan(3,j)+t2dz(j) enddo c energy calculation enenei=enenei+enei enddo edfanei=enenei print*, 'NRES: ',nres print*, 'EDFAN:',edfanei print*, 'COOR' do i=2,nres-1 write(*,'(a,i6,1x,3f12.5)')'DFANC:',i,c(1:3,i) enddo print*, 'NEI_FORCE' do i=2,nres-1 write(*,'(a,i6,1x,3f14.7)')'DFANG:',i,gdfan(1:3,i) enddo return end subroutine edfab(edfabeta) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.DFA' real*8 pai parameter(pai=3.14159265358979323846d0) C sheet variables real*8 bx(maxres),by(maxres),bz(maxres) real*8 vbet(maxres,maxres) real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres) real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12) real*8 vbeta,vbetp,vbetm real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, & c00,s00,ulnex,dnex real*8 dp45,dm45,w_beta common /sheca/ bx,by,bz common /shee/ vbeta,vbet,vbetp,vbetm common /shetf/ shetfx,shetfy,shetfz common /shef/ shefx, shefy, shefz common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, & c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta C End of sheet variables integer i,j double precision enebet enebet=0.0d0 bx=0.0d0;by=0.0d0;bz=0.0d0 shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 gdfab=0.0d0 c do i=1,nres do i=2,nres-1 bx(i-1)=c(1,i) by(i-1)=c(2,i) bz(i-1)=c(3,i) enddo dca=0.25d0**2 dshe=0.3d0**2 ulhb=5.0d0 uldhb=5.0d0 ulnex=cos(60.0d0/180.0d0*pai) dlhb=1.0d0 dldhb=1.0d0 dnex=0.3d0**2 c00=cos((1.0d0+10.0d0/180.0d0)*pai) s00=sin((1.0d0+10.0d0/180.0d0)*pai) w_beta=0.5d0 dp45=w_beta dm45=w_beta C END OF INITIALIZATION nca=nres-2 call angvectors(nca) call sheetforce(nca,wshet,dfaexp) c end of sheet energy and force c do j=1,nres do j=2,nres-1 shetfx(j-1)=shetfx(j-1)*beta_inc shetfy(j-1)=shetfy(j-1)*beta_inc shetfz(j-1)=shetfz(j-1)*beta_inc c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j) enddo vbeta=vbeta*beta_inc enebet=vbeta edfabeta=enebet c do j=1,nres do j=2,nres-1 gdfab(1,j)=gdfab(1,j)-shetfx(j-1) gdfab(2,j)=gdfab(2,j)-shetfy(j-1) gdfab(3,j)=gdfab(3,j)-shetfz(j-1) enddo print*, 'EDFAB:',edfabeta C print*, 'COOR' C do i=2,nres-1 C write(*,'(a,i6,1x,3f12.5)')'DFABC:',i,c(1:3,i) C enddo print*, 'BETA_FORCE' do i=2,nres-1 write(*,'(a,i6,1x,3f14.7)')'DFABG:',i,gdfab(1:3,i) enddo return end C------------------------------------------------------------------------------- subroutine angvectors(nca) c implicit real*4(a-h,o-z) implicit none integer nca integer maxca parameter(maxca=800) real*8 pai,zero parameter(PAI=3.14159265358979323846D0,zero=0.0d0) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 apx(maxca),apy(maxca),apz(maxca) real*8 apmx(maxca),apmy(maxca),apmz(maxca) real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca) real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca) real*8 atx(maxca),aty(maxca),atz(maxca) real*8 atmx(maxca),atmy(maxca),atmz(maxca) real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca) real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca) real*8 astx(maxca),asty(maxca),astz(maxca) real*8 astmx(maxca),astmy(maxca),astmz(maxca) real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca) real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca) real*8 sth(maxca) real*8 cph(maxca),cth(maxca) real*8 ulcos(maxca) real*8 p,c integer i, ip, ipp, ip3, j real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca) real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz real*8 cix, ciy, ciz, cipx, cipy, cipz real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz real*8 grpp,gx,gy,gz real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41 integer inb,nmax,iselect common /sheca/ bx,by,bz common /difvec/ rx, ry, rz common /ulang/ ulcos common /phys1/ inb,nmax,iselect common /phys4/ p,c common /kyori2/ dis common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy, & apmmz,apm3x,apm3y,apm3z common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy, & atmmz,atm3x,atm3y,atm3z common /coscos/ cph,cth common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy, & astmmz,astm3x,astm3y,astm3z common /sinsin/ sth C------------------------------------------------------------------------------- c write(*,*) 'inside angvectors' C initialize p=0.1d0 c=1.0d0 inb=nca cph=zero; cth=zero; sth=zero apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero astm3z=zero C end of initialize C r[x,y,z] calc and distance calculation rx=zero;ry=zero;rz=zero do i=1,inb do j=1,inb rx(i,j)=bx(j)-bx(i) ry(i,j)=by(j)-by(i) rz(i,j)=bz(j)-bz(i) dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2) c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i) c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i) c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i) c write(*,*) 'dis(i,j):',i,j,dis(i,j) enddo enddo c end of r[x,y,z] calc C cos calc do i=1,inb-2 ip=i+1 ipp=i+2 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp) $ +rz(i,ip)*rz(ip,ipp) ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp)) endif enddo c end of virtual bond angle c write(*,*) 'inside angvectors1' do i=1,inb-3 ip=i+1 ipp=i+2 ip3=i+3 rix=bx(ip)-bx(i) riy=by(ip)-by(i) riz=bz(ip)-bz(i) ripx=bx(ipp)-bx(ip) ripy=by(ipp)-by(ip) ripz=bz(ipp)-bz(ip) rippx=bx(ip3)-bx(ipp) rippy=by(ip3)-by(ipp) rippz=bz(ip3)-bz(ipp) gx=riy*ripz-riz*ripy gy=riz*ripx-rix*ripz gz=rix*ripy-riy*ripx gpx=ripy*rippz-ripz*rippy gpy=ripz*rippx-ripx*rippz gpz=ripx*rippy-ripy*rippx gpcrp_x=gpy*ripz-gpz*ripy gpcrp_y=gpz*ripx-gpx*ripz gpcrp_z=gpx*ripy-gpy*ripx d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2) gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy if(i.ge.2) then rimx=bx(i)-bx(i-1) rimy=by(i)-by(i-1) rimz=bz(i)-bz(i-1) gmx=rimy*riz-rimz*riy gmy=rimz*rix-rimx*riz gmz=rimx*riy-rimy*rix dgm=sqrt(gmx**2+gmy**2+gmz**2) dgm3=dgm**3 ggm=gmx*gx+gmy*gy+gmz*gz gmrp=gmx*ripx+gmy*ripy+gmz*ripz drim=dis(i-1,i) drim3=drim**3 gcr_x=gy*riz-gz*riy gcr_y=gz*rix-gx*riz gcr_z=gx*riy-gy*rix d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2) d_gcr3=d_gcr**3 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy endif c write(*,*) 'inside angvectors2' if(i.ge.3) then rimmx=bx(i-1)-bx(i-2) rimmy=by(i-1)-by(i-2) rimmz=bz(i-1)-bz(i-2) drimm=dis(i-2,i-1) gmmx=rimmy*rimz-rimmz*rimy gmmy=rimmz*rimx-rimmx*rimz gmmz=rimmx*rimy-rimmy*rimx dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2) dgmm3=dgmm**3 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz gmmr=gmmx*rix+gmmy*riy+gmmz*riz gmcrim_x=gmy*rimz-gmz*rimy gmcrim_y=gmz*rimx-gmx*rimz gmcrim_z=gmx*rimy-gmy*rimx d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2) d_gmcrim3=d_gmcrim**3 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy endif if(i.ge.4) then rim3x=bx(i-2)-bx(i-3) rim3y=by(i-2)-by(i-3) rim3z=bz(i-2)-bz(i-3) g3x=rim3y*rimmz-rim3z*rimmy g3y=rim3z*rimmx-rim3x*rimmz g3z=rim3x*rimmy-rim3y*rimmx dg30=sqrt(g3x**2+g3y**2+g3z**2) g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz g3rim_=g3x*rimx+g3y*rimy+g3z*rimz cc********************************************************************** gmmcrimm_x=gmmy*rimmz-gmmz*rimmy gmmcrimm_y=gmmz*rimmx-gmmx*rimmz gmmcrimm_z=gmmx*rimmy-gmmy*rimmx d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2) d_gmmcrimm3=d_gmmcrimm**3 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y endif dri=dis(i,i+1) drip=dis(i+1,i+2) dripp=dis(i+2,i+3) dri3=dri**3 dg=sqrt(gx**2+gy**2+gz**2) dgp=sqrt(gpx**2+gpy**2+gpz**2) dg3=dg**3 ggp=gx*gpx+gy*gpy+gz*gpz grpp=gx*rippx+gy*rippy+gz*rippz if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0 & .and.d_gpcrp.gt.0.0D0) then cph(i)=grpp/dg/dripp cth(i)=ggp/dg/dgp sth(i)=gpcrp__g/d_gpcrp/dg else c cph(i)=1.0D0 cth(i)=1.0D0 sth(i)=0.0D0 endif c write(*,*) 'inside angvectors3' if(dgp.gt.0.0D0.and.dg3.gt.0.0D0 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then d10=1.0D0/(dg*dgp) d11=ggp/(dg3*dgp) d12=1.0D0/(dg*dripp) d13=grpp/(dg3*dripp) sd10=1.0D0/(d_gpcrp*dg) sd11=gpcrp__g/(d_gpcrp*dg3) else d10=0.0D0 d11=0.0D0 d12=0.0D0 d13=0.0D0 sd10=0.0D0 sd11=0.0D0 endif atx(i)=(ripz*gpy-ripy*gpz)*d10 & -(gy*ripz-gz*ripy)*d11 aty(i)=(ripx*gpz-ripz*gpx)*d10 & -(gz*ripx-gx*ripz)*d11 atz(i)=(ripy*gpx-ripx*gpy)*d10 & -(gx*ripy-gy*ripx)*d11 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz & +ripy*gpy*ripx-gpx*ripz**2) & -sd11*(gy*ripz-gz*ripy) asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx & -gpy*ripx**2+gpz*ripy*ripz) & -sd11*(-gx*ripz+gz*ripx) astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2 & -gpz*ripy**2+ripz*gpx*ripx) & -sd11*(gx*ripy-gy*ripx) apx(i)=(ripz*rippy-ripy*rippz)*d12 & -(gy*ripz-gz*ripy)*d13 apy(i)=(ripx*rippz-ripz*rippx)*d12 & -(gz*ripx-gx*ripz)*d13 apz(i)=(ripy*rippx-ripx*rippy)*d12 & -(gx*ripy-gy*ripx)*d13 if(i.ge.2) then cix=bx(ip)-bx(i-1) ciy=by(ip)-by(i-1) ciz=bz(ip)-bz(i-1) cipx=bx(ipp)-bx(i) cipy=by(ipp)-by(i) cipz=bz(ipp)-bz(i) ripx=bx(ipp)-bx(ip) ripy=by(ipp)-by(ip) ripz=bz(ipp)-bz(ip) if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0 & .and.d_gcr3.gt.0.0D0) then d20=1.0D0/(dg*dgm) d21=ggm/(dgm3*dg) d22=ggm/(dgm*dg3) d23=1.0D0/(dgm*drip) d24=gmrp/(dgm3*drip) sd20=1.0D0/(d_gcr*dgm) sd21=gcr__gm/(d_gcr3*dgm) sd22=gcr__gm/(d_gcr*dgm3) else d20=0.0D0 d21=0.0D0 d22=0.0D0 d23=0.0D0 d24=0.0D0 sd20=0.0D0 sd21=0.0D0 sd22=0.0D0 endif atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20 & -(ciy*gmz-ciz*gmy)*d21 & +(ripy*gz-ripz*gy)*d22 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20 & -(ciz*gmx-cix*gmz)*d21 & +(ripz*gx-ripx*gz)*d22 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20 & -(cix*gmy-ciy*gmx)*d21 & +(ripx*gy-ripy*gx)*d22 cc********************************************************************** astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix & +gmz*gy+ripy*riy*gmx+riz*gx*ciz) & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz) & +gcr_z*(-ripz*rix+gy)) & -sd22*(-gmy*ciz+gmz*ciy) astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz & +riz*ripz*gmy) & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz) & -gcr_z*(ripz*riy+gx)) & -sd22*(gmx*ciz-gmz*cix) astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy & -riz*gx*cix) & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx) & +gcr_z*(ripy*riy+ripx*rix)) & -sd22*(-gmx*ciy+gmy*cix) cc********************************************************************** apmx(i)=(ciy*ripz-ripy*ciz)*d23 & -(ciy*gmz-ciz*gmy)*d24 apmy(i)=(ciz*ripx-ripz*cix)*d23 & -(ciz*gmx-cix*gmz)*d24 apmz(i)=(cix*ripy-ripx*ciy)*d23 & -(cix*gmy-ciy*gmx)*d24 endif if(i.ge.3) then if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0 & .and.d_gmcrim3.gt.0.0D0) then d30=1.0D0/(dgm*dgmm) d31=gmmgm/(dgm3*dgmm) d32=gmmgm/(dgm*dgmm3) d33=1.0D0/(dgmm*dri) d34=gmmr/(dgmm3*dri) d35=gmmr/(dgmm*dri3) sd30=1.0D0/(d_gmcrim*dgmm) sd31=gmcrim__gmm/(d_gmcrim3*dgmm) sd32=gmcrim__gmm/(d_gmcrim*dgmm3) else d30=0.0D0 d31=0.0D0 d32=0.0D0 d33=0.0D0 d34=0.0D0 d35=0.0D0 sd30=0.0D0 sd31=0.0D0 sd32=0.0D0 endif c write(*,*) 'inside angvectors4' cc********************************************************************** atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30 & -(ciy*gmz-ciz*gmy)*d31 & -(gmmy*rimmz-gmmz*rimmy)*d32 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30 & -(ciz*gmx-cix*gmz)*d31 & -(gmmz*rimmx-gmmx*rimmz)*d32 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30 & -(cix*gmy-ciy*gmx)*d31 & -(gmmx*rimmy-gmmy*rimmx)*d32 cc********************************************************************** astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy & -ciy*rimy*gmmx-rimz*gmx*rimmz) & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy) & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy)) & -sd32*(gmmy*rimmz-rimmy*gmmz) astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx & +gmz*rimy*rimmz-rimz*ciz*gmmy) & -sd31*(gmcrim_x*(cix*rimy-gmz) & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx)) & -sd32*(-gmmx*rimmz+rimmx*gmmz) astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy & +rimz*ciy*gmmy+rimz*gmx*rimmx) & -sd31*(gmcrim_x*(cix*rimz+gmy) & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx)) & -sd32*(gmmx*rimmy-rimmx*gmmy) c********************************************************************** apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33 & -(gmmy*rimmz-gmmz*rimmy)*d34 & +rix*d35 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33 & -(gmmz*rimmx-gmmx*rimmz)*d34 & +riy*d35 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33 & -(gmmx*rimmy-gmmy*rimmx)*d34 & +riz*d35 endif if(i.ge.4) then if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0 & .and.drim3.gt.0.0D0 & .and.d_gmmcrimm3.gt.0.0D0) then d40=1.0D0/(dg30*dgmm) d41=g3gmm/(dg30*dgmm3) d42=1.0D0/(dg30*drim) d43=g3rim_/(dg30*drim3) sd40=1.0D0/(dg30*d_gmmcrimm) sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30) else d40=0.0D0 d41=0.0D0 d42=0.0D0 d43=0.0D0 sd40=0.0D0 sd41=0.0D0 endif atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40 & -(gmmy*rimmz-gmmz*rimmy)*d41 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40 & -(gmmz*rimmx-gmmx*rimmz)*d41 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40 & -(gmmx*rimmy-gmmy*rimmx)*d41 cc********************************************************************** astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y & -g3z*rimmz*rimmx+rimmy**2*g3x) & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2) & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx) astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y & -rimmx*rimmy*g3x+rimmz**2*g3y) & -sd41*(-gmmcrimm_x*rimmx*rimmy & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx) astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z & +g3z*rimmx**2-rimmz*rimmy*g3y) & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz & +gmmcrimm_z*(rimmy**2+rimmx**2)) c********************************************************************** apm3x(i)=g3x*d42-rimx*d43 apm3y(i)=g3y*d42-rimy*d43 apm3z(i)=g3z*d42-rimz*d43 endif enddo c******************************************************************************* c write(*,*) 'inside angvectors5' do i=inb-2,inb rimx=bx(i)-bx(i-1) rimy=by(i)-by(i-1) rimz=bz(i)-bz(i-1) rimmx=bx(i-1)-bx(i-2) rimmy=by(i-1)-by(i-2) rimmz=bz(i-1)-bz(i-2) rim3x=bx(i-2)-bx(i-3) rim3y=by(i-2)-by(i-3) rim3z=bz(i-2)-bz(i-3) gmmx=rimmy*rimz-rimmz*rimy gmmy=rimmz*rimx-rimmx*rimz gmmz=rimmx*rimy-rimmy*rimx g3x=rim3y*rimmz-rim3z*rimmy g3y=rim3z*rimmx-rim3x*rimmz g3z=rim3x*rimmy-rim3y*rimmx dg30=sqrt(g3x**2+g3y**2+g3z**2) g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2) dgmm3=dgmm**3 drim=dis(i-1,i) drimm=dis(i-2,i-1) drim3=drim**3 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz cc********************************************************************** gmmcrimm_x=gmmy*rimmz-gmmz*rimmy gmmcrimm_y=gmmz*rimmx-gmmx*rimmz gmmcrimm_z=gmmx*rimmy-gmmy*rimmx d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2) d_gmmcrimm3=d_gmmcrimm**3 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0 & .and.drim3.gt.0.0D0 & .and.d_gmmcrimm3.gt.0.0D0) then d40=1.0D0/(dg30*dgmm) d41=g3gmm/(dg30*dgmm3) d42=1.0D0/(dg30*drim) d43=g3rim_/(dg30*drim3) sd40=1.0D0/(dg30*d_gmmcrimm) sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30) else d40=0.0D0 d41=0.0D0 d42=0.0D0 d43=0.0D0 sd40=0.0D0 sd41=0.0D0 endif atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40 & -(gmmy*rimmz-gmmz*rimmy)*d41 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40 & -(gmmz*rimmx-gmmx*rimmz)*d41 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40 & -(gmmx*rimmy-gmmy*rimmx)*d41 cc********************************************************************** astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y & -g3z*rimmz*rimmx+rimmy**2*g3x) & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2) & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx) astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y & -rimmx*rimmy*g3x+rimmz**2*g3y) & -sd41*(-gmmcrimm_x*rimmx*rimmy & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx) astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z & +g3z*rimmx**2-rimmz*rimmy*g3y) & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz & +gmmcrimm_z*(rimmy**2+rimmx**2)) cc********************************************************************** apm3x(i)=g3x*d42-rimx*d43 apm3y(i)=g3y*d42-rimy*d43 apm3z(i)=g3z*d42-rimz*d43 if(i.le.inb-1) then ip=i+1 rix=bx(ip)-bx(i) riy=by(ip)-by(i) riz=bz(ip)-bz(i) cix=bx(ip)-bx(i-1) ciy=by(ip)-by(i-1) ciz=bz(ip)-bz(i-1) gmx=rimy*riz-rimz*riy gmy=rimz*rix-rimx*riz gmz=rimx*riy-rimy*rix dgm=sqrt(gmx**2+gmy**2+gmz**2) dgm3=dgm**3 dri=dis(i,i+1) dri3=dri**3 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz gmmr=gmmx*rix+gmmy*riy+gmmz*riz gmcrim_x=gmy*rimz-gmz*rimy gmcrim_y=gmz*rimx-gmx*rimz gmcrim_z=gmx*rimy-gmy*rimx d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2) d_gmcrim3=d_gmcrim**3 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy if(dgm3.gt.0.0D0.and. & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0 & .and.d_gmcrim3.gt.0.0D0) then d30=1.0D0/(dgm*dgmm) d31=gmmgm/(dgm3*dgmm) d32=gmmgm/(dgm*dgmm3) d33=1.0D0/(dgmm*dri) d34=gmmr/(dgmm3*dri) d35=gmmr/(dgmm*dri3) sd30=1.0D0/(d_gmcrim*dgmm) sd31=gmcrim__gmm/(d_gmcrim3*dgmm) sd32=gmcrim__gmm/(d_gmcrim*dgmm3) else d30=0.0D0 d31=0.0D0 d32=0.0D0 d33=0.0D0 d34=0.0D0 d35=0.0D0 sd30=0.0D0 sd31=0.0D0 sd32=0.0D0 endif cc********************************************************************** atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30 & -(ciy*gmz-ciz*gmy)*d31 & -(gmmy*rimmz-gmmz*rimmy)*d32 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30 & -(ciz*gmx-cix*gmz)*d31 & -(gmmz*rimmx-gmmx*rimmz)*d32 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30 & -(cix*gmy-ciy*gmx)*d31 & -(gmmx*rimmy-gmmy*rimmx)*d32 cc********************************************************************** astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy & -ciy*rimy*gmmx-rimz*gmx*rimmz) & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy) & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy)) & -sd32*(gmmy*rimmz-rimmy*gmmz) astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx & +gmz*rimy*rimmz-rimz*ciz*gmmy) & -sd31*(gmcrim_x*(cix*rimy-gmz) & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx)) & -sd32*(-gmmx*rimmz+rimmx*gmmz) astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy & +rimz*ciy*gmmy+rimz*gmx*rimmx) & -sd31*(gmcrim_x*(cix*rimz+gmy) & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx)) & -sd32*(gmmx*rimmy-rimmx*gmmy) cc********************************************************************** apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33 & -(gmmy*rimmz-gmmz*rimmy)*d34 & +rix*d35 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33 & -(gmmz*rimmx-gmmx*rimmz)*d34 & +riy*d35 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33 & -(gmmx*rimmy-gmmy*rimmx)*d34 & +riz*d35 endif c write(*,*) 'inside angvectors6' if(i.eq.inb-2) then ipp=i+2 ripx=bx(ipp)-bx(ip) ripy=by(ipp)-by(ip) ripz=bz(ipp)-bz(ip) cipx=bx(ipp)-bx(i) cipy=by(ipp)-by(i) cipz=bz(ipp)-bz(i) gx=riy*ripz-riz*ripy gy=riz*ripx-rix*ripz gz=rix*ripy-riy*ripx ggm=gmx*gx+gmy*gy+gmz*gz gmrp=gmx*ripx+gmy*ripy+gmz*ripz dg=sqrt(gx**2+gy**2+gz**2) dg3=dg**3 drip=dis(i+1,i+2) gcr_x=gy*riz-gz*riy gcr_y=gz*rix-gx*riz gcr_z=gx*riy-gy*rix d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2) d_gcr3=d_gcr**3 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy if(dgm3.gt.0.0D0.and. & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0 & ) then d20=1.0D0/(dg*dgm) d21=ggm/(dgm3*dg) d22=ggm/(dgm*dg3) d23=1.0D0/(dgm*drip) d24=gmrp/(dgm3*drip) sd20=1.0D0/(d_gcr*dgm) sd21=gcr__gm/(d_gcr3*dgm) sd22=gcr__gm/(d_gcr*dgm3) else d20=0.0D0 d21=0.0D0 d22=0.0D0 d23=0.0D0 d24=0.0D0 sd20=0.0D0 sd21=0.0D0 sd22=0.0D0 endif atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20 & -(ciy*gmz-ciz*gmy)*d21 & +(ripy*gz-ripz*gy)*d22 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20 & -(ciz*gmx-cix*gmz)*d21 & +(ripz*gx-ripx*gz)*d22 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20 & -(cix*gmy-ciy*gmx)*d21 & +(ripx*gy-ripy*gx)*d22 cc********************************************************************** astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix & +gmz*gy+ripy*riy*gmx+riz*gx*ciz) & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz) & +gcr_z*(-ripz*rix+gy)) & -sd22*(-gmy*ciz+gmz*ciy) astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz & +riz*ripz*gmy) & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz) & -gcr_z*(ripz*riy+gx)) & -sd22*(gmx*ciz-gmz*cix) astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy & -riz*gx*cix) & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx) & +gcr_z*(ripy*riy+ripx*rix)) & -sd22*(-gmx*ciy+gmy*cix) cc********************************************************************** c apmx(i)=(ciy*ripz-ripy*ciz)*d23 & -(ciy*gmz-ciz*gmy)*d24 apmy(i)=(ciz*ripx-ripz*cix)*d23 & -(ciz*gmx-cix*gmz)*d24 apmz(i)=(cix*ripy-ripx*ciy)*d23 & -(cix*gmy-ciy*gmx)*d24 endif enddo return end c END of angvectors c------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- subroutine sheetforce(nca,wshet,dfaexp) implicit none C JYLEE c this should be matched with dfa.fcm integer maxca parameter(maxca=800) cc********************************************************************** integer nca integer i,k integer inb,nmax,iselect real*8 dfaexp(15001) real*8 vbeta,vbetp,vbetm real*8 shefx(maxca,12) real*8 shefy(maxca,12),shefz(maxca,12) real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca) real*8 vbet(maxca,maxca) real*8 wshet(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) common /sheca/ bx,by,bz common /phys1/ inb,nmax,iselect common /shef/ shefx,shefy,shefz common /shee/ vbeta,vbet,vbetp,vbetm common /shetf/ shetfx,shetfy,shetfz inb=nca do i=1,inb shetfx(i)=0.0D0 shetfy(i)=0.0D0 shetfz(i)=0.0D0 enddo do k=1,12 do i=1,inb shefx(i,k)=0.0D0 shefy(i,k)=0.0D0 shefz(i,k)=0.0D0 enddo enddo call sheetene(nca,wshet,dfaexp) call sheetforce1 887 format(a,1x,i6,3x,f12.8) 888 format(a,1x,i4,1x,i4,3x,f12.8) 889 format(a,1x,i4,3x,f12.8) !write(2,*) 'coord : ' do i=1,inb !write(2,887) 'bx:',i,bx(i) !write(2,887) 'by:',i,by(i) !write(2,887) 'bz:',i,bz(i) enddo !write(2,*) 'After sheetforce1' do i=1,inb do k=1,12 !write(2,888) 'shefx :',i,k,shefx(i,k) !write(2,888) 'shefy :',i,k,shefy(i,k) !write(2,888) 'shefz :',i,k,shefz(i,k) enddo enddo call sheetforce5 !write(2,*) 'After sheetforce5' do i=1,inb do k=1,12 !write(2,888) 'shefx :',i,k,shefx(i,k) !write(2,888) 'shefy :',i,k,shefy(i,k) !write(2,888) 'shefz :',i,k,shefz(i,k) enddo enddo call sheetforce6 !write(2,*) 'After sheetforce6' do i=1,inb do k=1,12 !write(2,888) 'shefx :',i,k,shefx(i,k) !write(2,888) 'shefy :',i,k,shefy(i,k) !write(2,888) 'shefz :',i,k,shefz(i,k) enddo enddo call sheetforce11 !write(2,*) 'After sheetforce11' do i=1,inb do k=1,12 !write(2,888) 'shefx :',i,k,shefx(i,k) !write(2,888) 'shefy :',i,k,shefy(i,k) !write(2,888) 'shefz :',i,k,shefz(i,k) enddo enddo call sheetforce12 !write(2,*) 'After sheetforce12' do i=1,inb do k=1,12 !write(2,888) 'shefx :',i,k,shefx(i,k) !write(2,888) 'shefy :',i,k,shefy(i,k) !write(2,888) 'shefz :',i,k,shefz(i,k) enddo enddo do i=1,inb do k=1,12 shetfx(i)=shetfx(i)+shefx(i,k) shetfy(i)=shetfy(i)+shefy(i,k) shetfz(i)=shetfz(i)+shefz(i,k) enddo enddo !write(2,*) 'Beta Finished' do i=1,inb !write(2,889) 'shetfx : ',i,shetfx(i) !write(2,889) 'shetfy : ',i,shetfy(i) !write(2,889) 'shetfz : ',i,shetfz(i) enddo return end C end sheetforce c------------------------------------------------------------------------------- subroutine sheetene(nca,wshet,dfaexp) implicit none integer maxca parameter(maxca=800) cc****************************************************************************** real*8 dfaexp(15001) real*8 dtmp1, dtmp2, dtmp3 real*8 vbet(maxca,maxca) real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 pin1(maxca,maxca),pin2(maxca,maxca) real*8 pin3(maxca,maxca),pin4(maxca,maxca) real*8 pina1(maxca,maxca),pina2(maxca,maxca) real*8 pina3(maxca,maxca),pina4(maxca,maxca) real*8 cph(maxca),cth(maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 ulcos(maxca) cc********************************************************************** real*8 astx(maxca),asty(maxca),astz(maxca) real*8 astmx(maxca),astmy(maxca),astmz(maxca) real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca) real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca) real*8 sth(maxca) real*8 wshet(maxca,maxca) real*8 dp45, dm45, w_beta real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb integer nca integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect real*8 uum, uup real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2 common /sheca/ bx,by,bz common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /coscos/ cph,cth common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, & c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4 common /shee/ vbeta,vbet,vbetp,vbetm common /ulang/ ulcos cc********************************************************************** common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy, & astmmz,astm3x,astm3y,astm3z common /sinsin/ sth real*8 r_pair_mat(maxca,maxca) common /beta_p/ r_pair_mat C------------------------------------------------------------------------------- r_pair_mat = 0.0d0 do i=1,inb do j=1,inb r_pair_mat(i,j)=wshet(i,j) c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j) enddo enddo c stop c vbeta=0.0D0 vbetp=0.0D0 vbetm=0.0D0 do i=1,inb-7 do j=i+4,inb-3 ip=i+1 ipp=i+2 jp=j+1 jpp=j+2 cc********************************************************************** y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2 y1=-0.5d0*y1/dca y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2 y2=-0.5d0*y2/dnex y=y1+y2 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp) $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp)) pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp) $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp)) pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp) $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp)) pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp) $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp)) yshe1=pin1(i,j)**2+pin2(i,j)**2 yshe1=-0.5d0*yshe1/dshe yshe2=pin3(i,j)**2+pin4(i,j)**2 yshe2=-0.5d0*yshe2/dshe C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i) C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i) C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i) C write(*,*) 'dis(i,j):',i,j,dis(i,j) C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp) C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp) C write(*,*) 'pin1:',pin1(i,j) C write(*,*) 'pin2:',pin2(i,j) C write(*,*) 'pin3:',pin3(i,j) C write(*,*) 'pin4:',pin4(i,j) C write(*,*) 'y:',y C write(*,*) 'yy1:',yy1 C write(*,*) 'yy2:',yy2 C write(*,*) 'yshe1:',yshe1 C write(*,*) 'yshe2:',yshe2 c dtmp1 = y+yy1+yshe1 dtmp2 = y+yy2+yshe2 dtmp3 = y+yy1+yy2+yshe1+yshe2 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3 C write(*,*)'2', y,yy1,yy2 C write(*,*)'3', yshe1,yshe2 if (dtmp3.le.-15.0d0) then c vbetap(i,j)=-dp45*exp(dtmp3) vbetap(i,j)=0.0d0 else vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1) end if if (dtmp1.le.-15.0d0) then c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1) vbetap1(i,j)=0.0d0 else vbetap1(i,j)=-r_pair_mat(i+1,j+1) $ *dfaexp(idint(-dtmp1*1000)+1) end if if (dtmp2.le.-15.0d0) then C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2) vbetap2(i,j)=0.0d0 else vbetap2(i,j)=-r_pair_mat(i+2,j+2) $ *dfaexp(idint(-dtmp2*1000)+1) end if c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2) c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1) c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2) ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1) ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2) yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp) $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp)) pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp) $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp)) pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp) $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp)) pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp) $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp)) yshe1=pina1(i,j)**2+pina2(i,j)**2 yshe1=-0.5d0*yshe1/dshe yshe2=pina3(i,j)**2+pina4(i,j)**2 yshe2=-0.5d0*yshe2/dshe C write(*,*) 'pina1:',pina1(i,j) C write(*,*) 'pina2:',pina2(i,j) C write(*,*) 'pina3:',pina3(i,j) C write(*,*) 'pina4:',pina4(i,j) C write(*,*) 'yshe1:',yshe1 C write(*,*) 'yshe2:',yshe2 C write(*,*) 'dshe:',dshe dtmp3=y+yy1+yy2+yshe1+yshe2 dtmp1=y+yy1+yshe1 dtmp2=y+yy2+yshe2 if(dtmp3 .le. -15.0d0) then c vbetam(i,j)=-dm45*exp(dtmp3) vbetam(i,j)=0.0d0 else vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1) end if if(dtmp1 .le. -15.0d0) then c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1) vbetam1(i,j)=0.0d0 else vbetam1(i,j)=-r_pair_mat(i+1,j+2) $ *dfaexp(idint(-dtmp1*1000)+1) end if if(dtmp2.le.-15.0d0) then c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2) vbetam2(i,j)=0.0d0 else vbetam2(i,j)=-r_pair_mat(i+2,j+1) $ *dfaexp(idint(-dtmp2*1000)+1) end if c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2) c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1) c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2) ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2) ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1) uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j) uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j) c write(*,*) 'uup,uum:', uup, uum c uup=vbetap1(i,j)+vbetap2(i,j) c uum=vbetam1(i,j)+vbetam2(i,j) vbet(i,j)=uup+uum vbetp=vbetp+uup vbetm=vbetm+uum vbeta=vbeta+vbet(i,j) c write(*,*) 'uup,uum:',uup,uum c write(*,*) 'vbetap(i,j):',vbetap(i,j) c write(*,*) 'vbetap1(i,j):',vbetap1(i,j) c write(*,*) 'vbetap2(i,j):',vbetap2(i,j) c write(*,*) 'vbetam(i,j):',vbetam(i,j) c write(*,*) 'vbetam1(i,j):',vbetam1(i,j) c write(*,*) 'vbetam2(i,j):',vbetam2(i,j) c write(*,*) 'uup:',uup c write(*,*) 'uum:',uum c write(*,*) 'vbetp:',vbetp c write(*,*) 'vbetm:',vbetm c write(*,*) 'vbet(i,j):',vbet(i,j) c stop enddo enddo ! do i=1,inb-7 ! do j=i+4,inb-3 ! write(*,*) 'I,J:', i,j ! write(*,*) 'vbetap(i,j):',vbetap(i,j) ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j) ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j) ! write(*,*) 'vbetam(i,j):',vbetam(i,j) ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j) ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j) ! write(*,*) 'vbet(i,j):',vbet(i,j) ! enddo ! enddo return end c------------------------------------------------------------------------------- subroutine sheetforce1 implicit none integer maxca parameter(maxca=800) cc********************************************************************** real*8 vbet(maxca,maxca) real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 cph(maxca),cth(maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 shefx(maxca,12) real*8 shefy(maxca,12),shefz(maxca,12) real*8 atx(maxca),aty(maxca),atz(maxca) real*8 atmx(maxca),atmy(maxca),atmz(maxca) real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca) real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca) real*8 apx(maxca),apy(maxca),apz(maxca) real*8 apmx(maxca),apmy(maxca),apmz(maxca) real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca) real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca) real*8 ulcos(maxca) real*8 astx(maxca),asty(maxca),astz(maxca) real*8 astmx(maxca),astmy(maxca),astmz(maxca) real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca) real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca) real*8 sth(maxca) real*8 w_beta,dp45, dm45 real*8 vbeta, vbetp, vbetm real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex integer inb,nmax,iselect common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /coscos/ cph,cth common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy, $ atmmz,atm3x,atm3y,atm3z common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy, $ apmmz,apm3x,apm3y,apm3z common /shef/ shefx,shefy,shefz common /shee/ vbeta,vbet,vbetp,vbetm common /ulang/ ulcos c c********************************************************************** common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy, $ astmmz,astm3x,astm3y,astm3z common /sinsin/ sth C-------------------------------------------------------------------------------- c local variables integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2 C-------------------------------------------------------------------------------- do i=4,inb-4 im3=i-3 imm=i-2 im=i-1 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca v1=0.0D0 do j=i+1,inb-3 v1=v1+vbet(im3,j) enddo cc1=(ulcos(imm)-ulnex)/dnex dmm=cc1/(dis(imm,im)*dis(im,i)) dmm__=cc1*ulcos(imm)/dis(im,i)**2 fx=rx(imm,im)*dmm-rx(im,i)*dmm__ fy=ry(imm,im)*dmm-ry(im,i)*dmm__ fz=rz(imm,im)*dmm-rz(im,i)*dmm__ fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1 shefx(i,1)=fx*v1 shefy(i,1)=fy*v1 shefz(i,1)=fz*v1 enddo do i=3,inb-5 imm=i-2 im=i-1 ip=i+1 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca v2=0.0D0 do j=i+2,inb-3 v2=v2+vbet(imm,j) enddo cc1=(ulcos(imm)-ulnex)/dnex cc2=(ulcos(im)-ulnex)/dnex dmm1=cc1/(dis(imm,im)*dis(im,i)) dmm2=cc2/(dis(im,i)*dis(i,ip)) dmm1__=cc1*ulcos(imm)/dis(im,i)**2 dmm2_1=cc2*ulcos(im)/dis(im,i)**2 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2 cc********************************************************************** fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2 shefx(i,2)=fx*v2 shefy(i,2)=fy*v2 shefz(i,2)=fz*v2 enddo do i=2,inb-6 im=i-1 ip=i+1 ipp=i+2 c3=(cth(im)*c00+sth(im)*s00-1)/dca v3=0.0D0 do j=i+3,inb-3 v3=v3+vbet(im,j) enddo cc2=(ulcos(im)-ulnex)/dnex cc3=(ulcos(i)-ulnex)/dnex dmm2=cc2/(dis(im,i)*dis(i,ip)) dmm3=cc3/(dis(i,ip)*dis(ip,ipp)) dmm2_1=cc2*ulcos(im)/dis(im,i)**2 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2 dmm3__=cc3*ulcos(i)/dis(i,ip)**2 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__ fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__ fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__ fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3 shefx(i,3)=fx*v3 shefy(i,3)=fy*v3 shefz(i,3)=fz*v3 enddo do i=1,inb-7 ip=i+1 ipp=i+2 c4=(cth(i)*c00+sth(i)*s00-1)/dca v4=0.0D0 do j=i+4,inb-3 v4=v4+vbet(i,j) enddo cc3=(ulcos(i)-ulnex)/dnex dmm3=cc3/(dis(i,ip)*dis(ip,ipp)) dmm3__=cc3*ulcos(i)/dis(i,ip)**2 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__ fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__ fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__ fx=fx+(atx(i)*c00+astx(i)*s00)*c4 fy=fy+(aty(i)*c00+asty(i)*s00)*c4 fz=fz+(atz(i)*c00+astz(i)*s00)*c4 shefx(i,4)=fx*v4 shefy(i,4)=fy*v4 shefz(i,4)=fz*v4 enddo do j=8,inb jm3=j-3 jmm=j-2 jm=j-1 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca v7=0.0D0 do i=1,j-7 v7=v7+vbet(i,jm3) enddo cc7=(ulcos(jmm)-ulnex)/dnex dmm=cc7/(dis(jmm,jm)*dis(jm,j)) dmm__=cc7*ulcos(jmm)/dis(jm,j)**2 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__ fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__ fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__ fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7 shefx(j,7)=fx*v7 shefy(j,7)=fy*v7 shefz(j,7)=fz*v7 enddo do j=7,inb-1 jm=j-1 jmm=j-2 jp=j+1 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca v8=0.0D0 do i=1,j-6 v8=v8+vbet(i,jmm) enddo cc7=(ulcos(jmm)-ulnex)/dnex cc8=(ulcos(jm)-ulnex)/dnex dmm7=cc7/(dis(jmm,jm)*dis(jm,j)) dmm8=cc8/(dis(jm,j)*dis(j,jp)) dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8 shefx(j,8)=fx*v8 shefy(j,8)=fy*v8 shefz(j,8)=fz*v8 enddo do j=6,inb-2 jm=j-1 jp=j+1 jpp=j+2 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca v9=0.0D0 do i=1,j-5 v9=v9+vbet(i,jm) enddo cc8=(ulcos(jm)-ulnex)/dnex cc9=(ulcos(j)-ulnex)/dnex dmm8=cc8/(dis(jm,j)*dis(j,jp)) dmm9=cc9/(dis(j,jp)*dis(jp,jpp)) dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2 dmm9__=cc9*ulcos(j)/dis(j,jp)**2 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__ fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__ fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__ fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9 shefx(j,9)=fx*v9 shefy(j,9)=fy*v9 shefz(j,9)=fz*v9 enddo do j=5,inb-3 jp=j+1 jpp=j+2 c10=(cth(j)*c00+sth(j)*s00-1)/dca v10=0.0D0 do i=1,j-4 v10=v10+vbet(i,j) enddo cc9=(ulcos(j)-ulnex)/dnex dmm9=cc9/(dis(j,jp)*dis(jp,jpp)) dmm9__=cc9*ulcos(j)/dis(j,jp)**2 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__ fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__ fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__ fx=fx+(atx(j)*c00+astx(j)*s00)*c10 fy=fy+(aty(j)*c00+asty(j)*s00)*c10 fz=fz+(atz(j)*c00+astz(j)*s00)*c10 shefx(j,10)=fx*v10 shefy(j,10)=fy*v10 shefz(j,10)=fz*v10 enddo return end c---------------------------------------------------------------------------- subroutine sheetforce5 implicit none integer maxca parameter(maxca=800) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 pin1(maxca,maxca),pin2(maxca,maxca) real*8 pin3(maxca,maxca),pin4(maxca,maxca) real*8 pina1(maxca,maxca),pina2(maxca,maxca) real*8 pina3(maxca,maxca),pina4(maxca,maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 shefx(maxca,12),shefy(maxca,12) real*8 shefz(maxca,12) real*8 dp45,dm45,w_beta real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex integer inb,nmax,iselect cc********************************************************************** common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4 common /shef/ shefx,shefy,shefz c******************************************************************************** c local variables integer i,imm,im,jp,jpp,j real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b c******************************************************************************** do i=3,inb-5 imm=i-2 im=i-1 do j=i+2,inb-3 jp=j+1 jpp=j+2 yy1=-(dis(i,jpp)-ulhb)/dlhb y1x=rx(jpp,i)/dis(i,jpp) y1y=ry(jpp,i)/dis(i,jpp) y1z=rz(jpp,i)/dis(i,jpp) y11x=yy1*y1x y11y=yy1*y1y y11z=yy1*y1z yy33=1.0D0/(dis(im,jp)*dis(im,i)) yyy3=pin1(imm,j)/(dis(im,i)**2) yy3=-pin1(imm,j)/dshe y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3 yy44=1.0D0/(dis(i,jpp)*dis(im,i)) yyy4a=pin3(imm,j)/(dis(i,jpp)**2) yyy4b=pin3(imm,j)/(dis(im,i)**2) yy4=-pin3(imm,j)/dshe y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp) $ -yyy4b*rx(im,i))*yy4 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp) $ -yyy4b*ry(im,i))*yy4 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp) $ -yyy4b*rz(im,i))*yy4 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp)) yyy5=pin4(imm,j)/(dis(i,jpp)**2) yy5=-pin4(imm,j)/dshe y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5 sx=y11x+y3x+y4x+y5x sy=y11y+y3y+y4y+y5y sz=y11z+y3z+y4z+y5z sx1=y3x sy1=y3y sz1=y3z sx2=y11x+y4x+y5x sy2=y11y+y4y+y5y sz2=y11z+y4z+y5z shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j) $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j) shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j) $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j) shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j) $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j) ! shefx(i,5)=shefx(i,5) ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j) ! shefy(i,5)=shefy(i,5) ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j) ! shefz(i,5)=shefz(i,5) ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j) yy6=-(dis(i,jp)-uldhb)/dldhb y6x=rx(jp,i)/dis(i,jp) y6y=ry(jp,i)/dis(i,jp) y6z=rz(jp,i)/dis(i,jp) y66x=yy6*y6x y66y=yy6*y6y y66z=yy6*y6z yy88=1.0D0/(dis(im,jpp)*dis(im,i)) yyy8=pina1(imm,j)/(dis(im,i)**2) yy8=-pina1(imm,j)/dshe y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8 yy99=1.0D0/(dis(jp,i)*dis(im,i)) yyy9a=pina3(imm,j)/(dis(jp,i)**2) yyy9b=pina3(imm,j)/(dis(im,i)**2) yy9=-pina3(imm,j)/dshe y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i) $ -yyy9b*rx(im,i))*yy9 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i) $ -yyy9b*ry(im,i))*yy9 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i) $ -yyy9b*rz(im,i))*yy9 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp)) yyy10=pina4(imm,j)/(dis(jp,i)**2) yy10=-pina4(imm,j)/dshe y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10 sx=y66x+y8x+y9x+y10x sy=y66y+y8y+y9y+y10y sz=y66z+y8z+y9z+y10z sx1=y8x sy1=y8y sz1=y8z sx2=y66x+y9x+y10x sy2=y66y+y9y+y10y sz2=y66z+y9z+y10z shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j) $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j) shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j) $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j) shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j) $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j) ! shefx(i,5)=shefx(i,5) ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j) ! shefy(i,5)=shefy(i,5) ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j) ! shefz(i,5)=shefz(i,5) ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j) enddo enddo return end c--------------------------------------------------------------------------c subroutine sheetforce6 implicit none integer maxca parameter(maxca=800) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 pin1(maxca,maxca),pin2(maxca,maxca) real*8 pin3(maxca,maxca),pin4(maxca,maxca) real*8 pina1(maxca,maxca),pina2(maxca,maxca) real*8 pina3(maxca,maxca),pina4(maxca,maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 shefx(maxca,12),shefy(maxca,12) real*8 shefz(maxca,12) real*8 dp45,dm45,w_beta real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex integer inb,nmax,iselect cc********************************************************************** common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4 common /shef/ shefx,shefy,shefz cc********************************************************************** C local variables integer i,imm,im,jp,jpp,j,ip real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b C******************************************************************************** do i=2,inb-6 ip=i+1 im=i-1 do j=i+3,inb-3 jp=j+1 jpp=j+2 yy1=-(dis(i,jp)-ulhb)/dlhb y1x=rx(jp,i)/dis(i,jp) y1y=ry(jp,i)/dis(i,jp) y1z=rz(jp,i)/dis(i,jp) y11x=yy1*y1x y11y=yy1*y1y y11z=yy1*y1z yy33=1.0D0/(dis(i,jp)*dis(i,ip)) yyy3a=pin1(im,j)/(dis(i,jp)**2) yyy3b=pin1(im,j)/(dis(i,ip)**2) yy3=-pin1(im,j)/dshe y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp) $ +yyy3b*rx(i,ip))*yy3 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp) $ +yyy3b*ry(i,ip))*yy3 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp) $ +yyy3b*rz(i,ip))*yy3 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp)) yyy4=pin2(im,j)/(dis(i,jp)**2) yy4=-pin2(im,j)/dshe y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip)) yyy5=pin3(im,j)/(dis(i,ip)**2) yy5=-pin3(im,j)/dshe y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5 sx=y11x+y3x+y4x+y5x sy=y11y+y3y+y4y+y5y sz=y11z+y3z+y4z+y5z sx1=y11x+y3x+y4x sy1=y11y+y3y+y4y sz1=y11z+y3z+y4z sx2=y5x sy2=y5y sz2=y5z shefx(i,6)=shefx(i,6)-sx*vbetap(im,j) $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j) shefy(i,6)=shefy(i,6)-sy*vbetap(im,j) $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j) shefz(i,6)=shefz(i,6)-sz*vbetap(im,j) $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j) ! shefx(i,6)=shefx(i,6) ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j) ! shefy(i,6)=shefy(i,6) ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j) ! shefz(i,6)=shefz(i,6) ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j) yy6=-(dis(jpp,i)-uldhb)/dldhb y6x=rx(jpp,i)/dis(jpp,i) y6y=ry(jpp,i)/dis(jpp,i) y6z=rz(jpp,i)/dis(jpp,i) y66x=yy6*y6x y66y=yy6*y6y y66z=yy6*y6z yy88=1.0D0/(dis(i,jpp)*dis(i,ip)) yyy8a=pina1(im,j)/(dis(i,jpp)**2) yyy8b=pina1(im,j)/(dis(i,ip)**2) yy8=-pina1(im,j)/dshe y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp) $ +yyy8b*rx(i,ip))*yy8 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp) $ +yyy8b*ry(i,ip))*yy8 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp) $ +yyy8b*rz(i,ip))*yy8 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp)) yyy9=pina2(im,j)/(dis(i,jpp)**2) yy9=-pina2(im,j)/dshe y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip)) yyy10=pina3(im,j)/(dis(i,ip)**2) yy10=-pina3(im,j)/dshe y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10 sx=y66x+y8x+y9x+y10x sy=y66y+y8y+y9y+y10y sz=y66z+y8z+y9z+y10z sx1=y66x+y8x+y9x sy1=y66y+y8y+y9y sz1=y66z+y8z+y9z sx2=y10x sy2=y10y sz2=y10z shefx(i,6)=shefx(i,6)-sx*vbetam(im,j) $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j) shefy(i,6)=shefy(i,6)-sy*vbetam(im,j) $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j) shefz(i,6)=shefz(i,6)-sz*vbetam(im,j) $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j) ! shefx(i,6)=shefx(i,6) ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j) ! shefy(i,6)=shefy(i,6) ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j) ! shefz(i,6)=shefz(i,6) ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j) enddo enddo return end c----------------------------------------------------------------------- subroutine sheetforce11 implicit none integer maxca parameter(maxca=800) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 pin1(maxca,maxca),pin2(maxca,maxca) real*8 pin3(maxca,maxca),pin4(maxca,maxca) real*8 pina1(maxca,maxca),pina2(maxca,maxca) real*8 pina3(maxca,maxca),pina4(maxca,maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 shefx(maxca,12),shefy(maxca,12) real*8 shefz(maxca,12) real*8 dp45,dm45,w_beta real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex integer inb,nmax,iselect cc********************************************************************** common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4 common /shef/ shefx,shefy,shefz C******************************************************************************** C local variables integer j,jm,jmm,ip,i,ipp real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8 C******************************************************************************** do j=7,inb-1 jm=j-1 jmm=j-2 do i=1,j-6 ip=i+1 ipp=i+2 yy1=-(dis(ipp,j)-ulhb)/dlhb y1x=rx(ipp,j)/dis(ipp,j) y1y=ry(ipp,j)/dis(ipp,j) y1z=rz(ipp,j)/dis(ipp,j) y11x=yy1*y1x y11y=yy1*y1y y11z=yy1*y1z yy33=1.0D0/(dis(ip,jm)*dis(jm,j)) yyy3=pin2(i,jmm)/(dis(jm,j)**2) yy3=-pin2(i,jmm)/dshe y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp)) yyy4=pin3(i,jmm)/(dis(ipp,j)**2) yy4=-pin3(i,jmm)/dshe y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4 yy55=1.0D0/(dis(ipp,j)*dis(jm,j)) yyy5a=pin4(i,jmm)/(dis(ipp,j)**2) yyy5b=pin4(i,jmm)/(dis(jm,j)**2) yy5=-pin4(i,jmm)/dshe y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j) $ -yyy5b*rx(jm,j))*yy5 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j) $ -yyy5b*ry(jm,j))*yy5 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j) $ -yyy5b*rz(jm,j))*yy5 sx=y11x+y3x+y4x+y5x sy=y11y+y3y+y4y+y5y sz=y11z+y3z+y4z+y5z sx1=y3x sy1=y3y sz1=y3z sx2=y11x+y4x+y5x sy2=y11y+y4y+y5y sz2=y11z+y4z+y5z shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm) $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm) shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm) $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm) shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm) $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm) ! shefx(j,11)=shefx(j,11) ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm) ! shefy(j,11)=shefy(j,11) ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm) ! shefz(j,11)=shefz(j,11) ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm) yy6=-(dis(ip,j)-uldhb)/dldhb y6x=rx(ip,j)/dis(ip,j) y6y=ry(ip,j)/dis(ip,j) y6z=rz(ip,j)/dis(ip,j) y66x=yy6*y6x y66y=yy6*y6y y66z=yy6*y6z yy88=1.0D0/(dis(ip,j)*dis(ip,ipp)) yyy8=pina1(i,jmm)/(dis(ip,j)**2) yy8=-pina1(i,jmm)/dshe y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8 yy99=1.0D0/(dis(ip,j)*dis(jm,j)) yyy9a=pina2(i,jmm)/(dis(ip,j)**2) yyy9b=pina2(i,jmm)/(dis(jm,j)**2) yy9=-pina2(i,jmm)/dshe y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j) $ -yyy9b*rx(jm,j))*yy9 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j) $ -yyy9b*ry(jm,j))*yy9 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j) $ -yyy9b*rz(jm,j))*yy9 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j)) yyy10=pina4(i,jmm)/(dis(jm,j)**2) yy10=-pina4(i,jmm)/dshe y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10 sx=y66x+y8x+y9x+y10x sy=y66y+y8y+y9y+y10y sz=y66z+y8z+y9z+y10z sx1=y66x+y8x+y9x sy1=y66y+y8y+y9y sz1=y66z+y8z+y9z sx2=y10x sy2=y10y sz2=y10z shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm) $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm) shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm) $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm) shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm) $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm) ! shefx(j,11)=shefx(j,11) ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm) ! shefy(j,11)=shefy(j,11) ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm) ! shefz(j,11)=shefz(j,11) ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm) enddo enddo return end c----------------------------------------------------------------------- subroutine sheetforce12 implicit none integer maxca parameter(maxca=800) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca) real*8 pin1(maxca,maxca),pin2(maxca,maxca) real*8 pin3(maxca,maxca),pin4(maxca,maxca) real*8 pina1(maxca,maxca),pina2(maxca,maxca) real*8 pina3(maxca,maxca),pina4(maxca,maxca) real*8 rx(maxca,maxca) real*8 ry(maxca,maxca),rz(maxca,maxca) real*8 bx(maxca),by(maxca),bz(maxca) real*8 dis(maxca,maxca) real*8 shefx(maxca,12),shefy(maxca,12) real*8 shefz(maxca,12) real*8 dp45,dm45,w_beta real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex integer inb,nmax,iselect cc********************************************************************** common /phys1/ inb,nmax,iselect common /kyori2/ dis common /difvec/ rx,ry,rz common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb, $ c00,s00,ulnex,dnex common /sheconst/ dp45,dm45,w_beta common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4 common /shef/ shefx,shefy,shefz cc********************************************************************** C local variables integer j,jm,jmm,ip,i,ipp,jp real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8 !c*************************************************************************c do j=6,inb-2 jp=j+1 jm=j-1 do i=1,j-5 ip=i+1 ipp=i+2 yy1=-(dis(ip,j)-ulhb)/dlhb y1x=rx(ip,j)/dis(ip,j) y1y=ry(ip,j)/dis(ip,j) y1z=rz(ip,j)/dis(ip,j) y11x=y1x*yy1 y11y=y1y*yy1 y11z=y1z*yy1 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp)) yyy3=pin1(i,jm)/(dis(ip,j)**2) yy3=-pin1(i,jm)/dshe y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3 yy44=1.0D0/(dis(ip,j)*dis(j,jp)) yyy4a=pin2(i,jm)/(dis(ip,j)**2) yyy4b=pin2(i,jm)/(dis(j,jp)**2) yy4=-pin2(i,jm)/dshe y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j) $ +yyy4b*rx(j,jp))*yy4 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j) $ +yyy4b*ry(j,jp))*yy4 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j) $ +yyy4b*rz(j,jp))*yy4 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp)) yyy5=pin4(i,jm)/(dis(j,jp)**2) yy5=-pin4(i,jm)/dshe y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5 sx=y11x+y3x+y4x+y5x sy=y11y+y3y+y4y+y5y sz=y11z+y3z+y4z+y5z sx1=y11x+y3x+y4x sy1=y11y+y3y+y4y sz1=y11z+y3z+y4z sx2=y5x sy2=y5y sz2=y5z shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm) $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm) shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm) $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm) shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm) $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm) ! shefx(j,12)=shefx(j,12) ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm) ! shefy(j,12)=shefy(j,12) ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm) ! shefz(j,12)=shefz(j,12) ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm) yy6=-(dis(ipp,j)-uldhb)/dldhb y6x=rx(ipp,j)/dis(ipp,j) y6y=ry(ipp,j)/dis(ipp,j) y6z=rz(ipp,j)/dis(ipp,j) y66x=yy6*y6x y66y=yy6*y6y y66z=yy6*y6z yy88=1.0D0/(dis(ip,jp)*dis(j,jp)) yyy8=pina2(i,jm)/(dis(j,jp)**2) yy8=-pina2(i,jm)/dshe y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp)) yyy9=pina3(i,jm)/(dis(j,ipp)**2) yy9=-pina3(i,jm)/dshe y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp)) yyy10a=pina4(i,jm)/(dis(j,ipp)**2) yyy10b=pina4(i,jm)/(dis(j,jp)**2) yy10=-pina4(i,jm)/dshe y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp) $ +yyy10b*rx(j,jp))*yy10 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp) $ +yyy10b*ry(j,jp))*yy10 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp) $ +yyy10b*rz(j,jp))*yy10 sx=y66x+y8x+y9x+y10x sy=y66y+y8y+y9y+y10y sz=y66z+y8z+y9z+y10z sx1=y8x sy1=y8y sz1=y8z sx2=y66x+y9x+y10x sy2=y66y+y9y+y10y sz2=y66z+y9z+y10z shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm) $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm) shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm) $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm) shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm) $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm) ENDDO ENDDO RETURN END C===============================================================================