+++ /dev/null
- subroutine init_dfa_vars
-
- include 'DIMENSIONS'
- include 'COMMON.INTERACT'
- 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 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!
-c dfaexp = 0.0d0
-c do ii = 1, 15001
-c dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
-c end do
-
- ishiftca=nnt-1
- ilastca=nct
-
- print *,'ishiftca=',ishiftca,'ilastca=',ilastca
-
- 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, *) 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)
-
-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, *) 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, *) 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, *) 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,*) 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
-
- edfadis=0
- gdfad=0.0d0
-
- do i=1, idfadis
-
- iatm1=idislis(1,i)+ishiftca
- iatm2=idislis(2,i)+ishiftca
- 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
-
- do j=1,idisnum(i)
-
- dd = dist-fdist(i,j)
- dtmp = dd*dd/ckk
- if (dtmp.ge.15.0d0) then
- texp = 0.0d0
- else
-c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
- texp = exp(-dtmp)/sckk
- endif
-
- ep=ep+sccdist(i,j)*texp
- fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
- scc=scc+sccdist(i,j)
-C write(*,'(2i8,6f12.5)') i, j, dist,
-C & fdist(i,j), ep, fp, sccdist(i,j), scc
-
- enddo
-
- ep = -ep/scc
- fp = fp/scc
-
-
-c IF(ABS(EP).lt.1.0d-20)THEN
-c EP=0.0D0
-c ENDIF
-c IF (ABS(FP).lt.1.0d-20) THEN
-c FP=0.0D0
-c ENDIF
-
- edfadis=edfadis+ep*dis_inc*wwdist
-
- gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
- gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
- gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
-
- gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
- gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
- gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
-
- 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)+ishiftca
- enddo
-
-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
-c ps_tmp = dfaexp(idint(dtmp*1000)+1)
- ps_tmp = exp(-dtmp)
- 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)=ithelis(iii,i)+ishiftca
- enddo
-
-
-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
-c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
- th_tmp = exp(-dtmp)
- 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
-
- 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,dtmp
-
- 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)+ishiftca
-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
-
- do j = ishiftca+1, ilastca
-
- 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
-
- 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
-c ntmp = dfaexp( idint(dtmp*1000) + 1 )
- ntmp = exp(-dtmp)
- 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
-
-c if (abs(enei).lt.1.0d-20)then
-c enei=0.0d0
-c endif
-c if (abs(tmp_fnei).lt.1.0d-20) then
-c tmp_fnei=0.0d0
-c endif
-
-c force calculation
- t1dx=t1dx*tmp_fnei
- t1dy=t1dy*tmp_fnei
- t1dz=t1dz*tmp_fnei
-
- do j=ishiftca+1,ilastca
- 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=ishiftca+1,ilastca
- 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
-
- 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)
- parameter (maxca=800)
-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
-
- real*8 cph(maxca),cth(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 sth(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 atxnum(maxca),atynum(maxca),atznum(maxca),
- & astxnum(maxca),astynum(maxca),astznum(maxca),
- & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
- & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
- & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
- & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
- & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
- & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
- & cth_orig(maxca),sth_orig(maxca)
-
- 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
-
- common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
- $ atmmz,atm3x,atm3y,atm3z
- common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
- $ astmmz,astm3x,astm3y,astm3z
-
- common /coscos/ cph,cth
- common /sinsin/ sth
-
-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
-
- do i=ishiftca+1,ilastca
- bx(i-ishiftca)=c(1,i)
- by(i-ishiftca)=c(2,i)
- bz(i-ishiftca)=c(3,i)
- enddo
-
-c do i=1,ilastca-ishiftca
-c read(99,*) bx(i),by(i),bz(i)
-c enddo
-c close(99)
-
- 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=ilastca-ishiftca
-
- call angvectors(nca)
- call sheetforce(nca,wshet)
-
-c end of sheet energy and force
-
- do j=1,nca
- shetfx(j)=shetfx(j)*beta_inc
- shetfy(j)=shetfy(j)*beta_inc
- shetfz(j)=shetfz(j)*beta_inc
-c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
- enddo
-
- vbeta=vbeta*beta_inc
- enebet=vbeta
- edfabeta=enebet
-
- do j=1,nca
- gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
- gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
- gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
- enddo
-
-#ifdef DEBUG1
- do j=1,nca
- write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
- enddo
-
-
- gdfab=0
- dinc=0.001
- do j=1,nca
- cth_orig(j)=cth(j)
- sth_orig(j)=sth(j)
- enddo
-
- do j=1,nca
-
- bx(j)=bx(j)+dinc
- call angvectors(nca)
- bx(j)=bx(j)-2*dinc
- call angvectors(nca)
- atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
- astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
- if (j.gt.1) then
- atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
- astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
- endif
- if (j.gt.2) then
- atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
- astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
- endif
- if (j.gt.3) then
- atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
- astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
- endif
- bx(j)=bx(j)+dinc
- by(j)=by(j)+dinc
- call angvectors(nca)
- by(j)=by(j)-2*dinc
- call angvectors(nca)
- by(j)=by(j)+dinc
- atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
- astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
- if (j.gt.1) then
- atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
- astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
- endif
- if (j.gt.2) then
- atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
- astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
- endif
- if (j.gt.3) then
- atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
- astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
- endif
-
- bz(j)=bz(j)+dinc
- call angvectors(nca)
- bz(j)=bz(j)-2*dinc
- call angvectors(nca)
- bz(j)=bz(j)+dinc
-
- atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
- astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
- if (j.gt.1) then
- atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
- astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
- endif
- if (j.gt.2) then
- atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
- astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
- endif
- if (j.gt.3) then
- atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
- astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
- endif
-
- enddo
-
- do i=1,nca
- write (*,'(2i5,a2,6f10.5)')
- & i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
- & astxnum(i),astx(i),astxnum(i)/astx(i),
- & i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
- & astynum(i),asty(i),astynum(i)/asty(i),
- & i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
- & astznum(i),astz(i),astznum(i)/astz(i),
- & i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
- & astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
- & i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
- & astmynum(i),astmy(i),astmynum(i)/astmy(i),
- & i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
- & astmznum(i),astmz(i),astmznum(i)/astmz(i),
- & i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
- & astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
- & i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
- & astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
- & i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
- & astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
- & i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
- & astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
- & i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
- & astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
- & i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
- & astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
- & i,0," ",cth_orig(i),sth_orig(i)
- enddo
-
-
- gdfab=0
- dinc=0.001
-
- do j=1,nca
-
- bx(j)=bx(j)+dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta1=vbeta*beta_inc
- bx(j)=bx(j)-2*dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta2=vbeta*beta_inc
- gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
- bx(j)=bx(j)+dinc
-
- by(j)=by(j)+dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta1=vbeta*beta_inc
- by(j)=by(j)-2*dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta2=vbeta*beta_inc
- gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
- by(j)=by(j)+dinc
-
- bz(j)=bz(j)+dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta1=vbeta*beta_inc
- bz(j)=bz(j)-2*dinc
- call angvectors(nca)
- call sheetforce(nca,wshet)
- vbeta2=vbeta*beta_inc
- gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
- bz(j)=bz(j)+dinc
-
-
- enddo
-
-
- call angvectors(nca)
- call sheetforce(nca,wshet)
- do j=1,nca
- shetfx(j)=shetfx(j)*beta_inc
- shetfy(j)=shetfy(j)*beta_inc
- shetfz(j)=shetfz(j)*beta_inc
- enddo
-
-
- write(*,*) 'xyz analytical and numerical gradient'
- do j=1,nca
- write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
- & ,(-gdfab(i,j),i=1,3)
- enddo
-
- do j=1,nca
- write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
- & shetfy(j)/gdfab(2,j),
- & shetfz(j)/gdfab(3,j)
- enddo
-
- stop
-#endif
-
- 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'
-crc do i=1,inb-3
- do i=1,inb
- 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*rimmy)
-
-c & +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'
-
-c do i=inb-2,inb
- do i=1,0
- 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)
- 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
-
-c 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)
- 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)
- implicit none
- integer maxca
- parameter(maxca=800)
-cc******************************************************************************
-
-c 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)
-ci integer istrand(maxca,maxca)
-ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
-ci common /shetest/ istrand,istrand_p,istrand_m
- 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
-
-cdebug y2=0
-
- y=y1+y2
-
-ci if(y.ge.-4) then
-ci istrand(i,j)=1
-ci else
-ci istrand(i,j)=0
-ci endif
-
-ci if(istrand(i,j).eq.1) then
-
- 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
-
-ci if((yshe1+yshe2).ge.-4) then
-ci istrand_p(i,j)=1
-ci else
-ci istrand_p(i,j)=0
-ci endif
-
-
-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
-
-ci if (istrand_p(i,j).eq.1) then
-
-cd yy1=0
-cd yy2=0
-cd yshe1=0
-cd yshe2=0
- 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
-
-cc if (dtmp3.le.-35.0d0) then
-c vbetap(i,j)=-dp45*exp(dtmp3)
-cc vbetap(i,j)=0.0d0
-cc else
-c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
- vbetap(i,j)=-dp45*exp(dtmp3)
-cc end if
-
-cc if (dtmp1.le.-35.0d0) then
-c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
-cc vbetap1(i,j)=0.0d0
-cc else
-c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
-c $ *dfaexp(idint(-dtmp1*1000)+1)
- vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
-cc end if
-
-cc if (dtmp2.le.-35.0d0) then
-C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
-cc vbetap2(i,j)=0.0d0
-cc else
-c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
-c $ *dfaexp(idint(-dtmp2*1000)+1)
- vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
-cc 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)
-
-ci elseif (istrand_p(i,j).eq.0)then
-ci vbetap(i,j)=0
-ci vbetap1(i,j)=0
-ci vbetap2(i,j)=0
-ci endif
-
- 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
-
-ci if((yshe1+yshe2).ge.-4) then
-ci istrand_m(i,j)=1
-ci else
-ci istrand_m(i,j)=0
-ci endif
-
-
-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
-
-ci if (istrand_m(i,j).eq.1) then
-
-cd yy1=0
-cd yy2=0
-cd yshe1=0
-cd yshe2=0
-
- dtmp3=y+yy1+yy2+yshe1+yshe2
- dtmp1=y+yy1+yshe1
- dtmp2=y+yy2+yshe2
-
-cc if(dtmp3 .le. -35.0d0) then
-c vbetam(i,j)=-dm45*exp(dtmp3)
-cc vbetam(i,j)=0.0d0
-cc else
-c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
- vbetam(i,j)=-dm45*exp(dtmp3)
-cc end if
-
-cc if(dtmp1 .le. -35.0d0) then
-c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
-cc vbetam1(i,j)=0.0d0
-cc else
-c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
-c $ *dfaexp(idint(-dtmp1*1000)+1)
- vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
-cc end if
-
-cc if(dtmp2.le.-35.0d0) then
-c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
-cc vbetam2(i,j)=0.0d0
-cc else
-c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
-c $ *dfaexp(idint(-dtmp2*1000)+1)
- vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
-cc end if
-
-ci elseif (istrand_m(i,j).eq.0)then
-ci vbetam(i,j)=0
-ci vbetam1(i,j)=0
-ci vbetam2(i,j)=0
-ci endif
-
-
-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)
-
-ci elseif(istrand(i,j).eq.0)then
-ci vbet(i,j)=0
-ci endif
-
-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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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
-cd fx=0
-cd fy=0
-cd fz=0
- 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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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
-cd fx=0
-cd fy=0
-cd fz=0
- 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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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__
-cd fx=0
-cd fy=0
-cd fz=0
- 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
-ci integer istrand(maxca,maxca)
-ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
-ci common /shetest/ istrand,istrand_p,istrand_m
-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
-
-ci if(istrand(imm,j).eq.1
-ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
-
-
- 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)
-
-ci endif
-
- 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
-ci integer istrand(maxca,maxca)
-ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
-ci common /shetest/ istrand,istrand_p,istrand_m
-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
-
-ci if(istrand(im,j).eq.1
-ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
-
-
- 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)
-
-ci endif
-
- 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
-ci integer istrand(maxca,maxca)
-ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
-ci common /shetest/ istrand,istrand_p,istrand_m
-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
-
-ci if(istrand(i,jmm).eq.1
-ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
-
-
- 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)
-
-ci endif
-
- 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
-ci integer istrand(maxca,maxca)
-ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
-ci common /shetest/ istrand,istrand_p,istrand_m
-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
-
-ci if(istrand(i,jm).eq.1
-ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
-
-
- 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)
-
-ci endif
-
- ENDDO
- ENDDO
-
- RETURN
- END
-C===============================================================================
projects / unres.git / blobdiff