source/cluster/wham/src/dfa.F

   1       subroutine init_dfa_vars
   2
   3       include 'DIMENSIONS'
   4       include 'COMMON.INTERACT'
   5       include 'COMMON.DFA'
   6
   7       integer ii
   8
   9 C     Number of restraints
  10       idisnum = 0
  11       iphinum = 0
  12       ithenum = 0
  13       ineinum = 0
  14
  15       idislis = 0
  16       iphilis = 0
  17       ithelis = 0
  18       ineilis = 0
  19       jneilis = 0
  20       jneinum = 0
  21       kshell  = 0
  22       fnei    = 0
  23 C     For beta
  24       nca     = 0
  25       icaidx  = 0
  26
  27 C     real variables
  28 CC    WEIGHTS for each min
  29       sccdist = 0.0d0
  30       fdist   = 0.0d0
  31       sccphi  = 0.0d0
  32       sccthe  = 0.0d0
  33       sccnei  = 0.0d0
  34       fphi1   = 0.0d0
  35       fphi2   = 0.0d0
  36       fthe1   = 0.0d0
  37       fthe2   = 0.0d0
  38 C     energies
  39       edfatot = 0.0d0
  40       edfadis = 0.0d0
  41       edfaphi = 0.0d0
  42       edfathe = 0.0d0
  43       edfanei = 0.0d0
  44       edfabet = 0.0d0
  45 C     weights for each E term
  46 C     these should be identical with
  47       dis_inc = 0.0d0
  48       phi_inc = 0.0d0
  49       the_inc = 0.0d0
  50       nei_inc = 0.0d0
  51       beta_inc = 0.0d0
  52       wshet   = 0.0d0
  53 C     precalculate exp table!
  54 c      dfaexp  = 0.0d0
  55 c      do ii = 1, 15001
  56 c         dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
  57 c      end do
  58
  59       ishiftca=nnt-1
  60       ilastca=nct
  61
  62       print *,'ishiftca=',ishiftca,'ilastca=',ilastca
  63
  64       return
  65       end
  66
  67
  68       subroutine read_dfa_info
  69 C
  70 C     read fragment informations
  71 C
  72       implicit real*8 (a-h,o-z)
  73       include 'DIMENSIONS'
  74       include 'COMMON.IOUNITS'
  75       include 'COMMON.CHAIN'
  76       include 'COMMON.DFA'
  77
  78
  79 C     NOTE THAT FILENAMES are FIXED, CURRENTLY!!
  80 C     THIS SHOULD BE MODIFIED!!
  81
  82       character*320 buffer
  83       integer iodfa
  84       parameter(iodfa=89)
  85
  86       integer i, j, nval
  87       integer ica1, ica2,ica3,ica4,ica5
  88       integer ishell, inca, itmp,iitmp
  89       double precision wtmp
  90 C
  91 C     READ DISTANCE
  92 C
  93       open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
  94       goto 34
  95  33   write(iout,'(a)') 'Error opening dist_dfa.dat file'
  96       stop
  97  34   continue
  98       write(iout,'(a)') 'dist_dfa.dat is opened!'
  99 C     read title
 100       read(iodfa, '(a)') buffer
 101 C     read number of restraints
 102       read(iodfa, *) IDFADIS
 103       read(iodfa, *) dis_inc
 104       do i=1, idfadis
 105          read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
 106
 107          idisnum(i)=nval
 108          idislis(1,i)=ica1
 109          idislis(2,i)=ica2
 110
 111          do j=1, nval
 112             read(iodfa,*) tmp
 113             fdist(i,j) = tmp
 114          enddo
 115
 116          do j=1, nval
 117             read(iodfa,*) tmp
 118             sccdist(i,j) = tmp
 119          enddo
 120
 121       enddo
 122       close(iodfa)
 123
 124 C     READ ANGLE RESTRAINTS
 125 C     PHI RESTRAINTS
 126       open(iodfa, file='phi_dfa.dat',status='old',err=35)
 127       goto 36
 128  35   write(iout,'(a)') 'Error opening dist_dfa.dat file'
 129       stop
 130
 131  36   continue
 132       write(iout,'(a)') 'phi_dfa.dat is opened!'
 133
 134 C     READ TITLE
 135       read(iodfa, '(a)') buffer
 136 C     READ NUMBER OF RESTRAINTS
 137       READ(iodfa, *) IDFAPHI
 138       read(iodfa,*) phi_inc
 139       do i=1, idfaphi
 140          read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
 141
 142          iphinum(i)=nval
 143
 144          iphilis(1,i)=ica1
 145          iphilis(2,i)=ica2
 146          iphilis(3,i)=ica3
 147          iphilis(4,i)=ica4
 148          iphilis(5,i)=ica5
 149
 150          do j=1, nval
 151             read(iodfa,*) tmp1,tmp2
 152             fphi1(i,j) = tmp1
 153             fphi2(i,j) = tmp2
 154          enddo
 155
 156          do j=1, nval
 157             read(iodfa,*) tmp
 158             sccphi(i,j) = tmp
 159          enddo
 160
 161       enddo
 162       close(iodfa)
 163
 164 C     THETA RESTRAINTS
 165       open(iodfa, file='theta_dfa.dat',status='old',err=41)
 166       goto 42
 167  41   write(iout,'(a)') 'Error opening dist_dfa.dat file'
 168       stop
 169  42   continue
 170       write(iout,'(a)') 'theta_dfa.dat is opened!'
 171 C     READ TITLE
 172       read(iodfa, '(a)') buffer
 173 C     READ NUMBER OF RESTRAINTS
 174       READ(iodfa, *) IDFATHE
 175       read(iodfa,*) the_inc
 176
 177       do i=1, idfathe
 178          read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
 179
 180          ithenum(i)=nval
 181
 182          ithelis(1,i)=ica1
 183          ithelis(2,i)=ica2
 184          ithelis(3,i)=ica3
 185          ithelis(4,i)=ica4
 186          ithelis(5,i)=ica5
 187
 188          do j=1, nval
 189             read(iodfa,*) tmp1,tmp2
 190             fthe1(i,j) = tmp1
 191             fthe2(i,j) = tmp2
 192          enddo
 193
 194          do j=1, nval
 195             read(iodfa,*) tmp
 196             sccthe(i,j) = tmp
 197          enddo
 198
 199       enddo
 200       close(iodfa)
 201 C     END of READING ANGLE RESTRAINT!
 202
 203 C     NUMBER OF NEIGHBOR CAs
 204       open(iodfa,file='nei_dfa.dat',status='old',err=37)
 205       goto 38
 206  37   write(iout,'(a)') 'Error opening nei_dfa.dat file'
 207       stop
 208  38   continue
 209       write(iout,'(a)') 'nei_dfa.dat is opened!'
 210 C     READ TITLE
 211       read(iodfa, '(a)') buffer
 212 C     READ NUMBER OF RESTRAINTS
 213       READ(iodfa, *) idfanei
 214       read(iodfa,*) nei_inc
 215
 216       do i=1, idfanei
 217          read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
 218
 219          ineilis(i)=ica1
 220          kshell(i)=ishell
 221          ineinum(i)=nval
 222
 223          do j=1, nval
 224             read(iodfa,*) inca
 225             fnei(i,j) = inca
 226 C            write(*,*) 'READ NEI:',i,j,fnei(i,j)
 227          enddo
 228
 229          do j=1, nval
 230             read(iodfa,*) tmp
 231             sccnei(i,j) = tmp
 232          enddo
 233
 234       enddo
 235       close(iodfa)
 236 C     END OF NEIGHBORING CA
 237
 238 C     READ BETA RESTRAINT
 239       open(iodfa, file='beta_dfa.dat',status='old',err=39)
 240       goto 40
 241  39   write(iout,'(a)') 'Error opening beta_dfa.dat file'
 242       stop
 243  40   continue
 244       write(iout,'(a)') 'beta_dfa.dat is opened!'
 245
 246       read(iodfa,'(a)') buffer
 247       read(iodfa,*) itmp
 248       read(iodfa,*) beta_inc
 249
 250       do i=1,itmp
 251          read(iodfa,*) ica1, iitmp
 252          do j=1,itmp
 253             read(iodfa,*) wtmp
 254             wshet(i,j) =  wtmp
 255 c            write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
 256          enddo
 257       enddo
 258
 259       close(iodfa)
 260 C     END OF BETA RESTRAINT
 261
 262       return
 263       END
 264
 265       subroutine edfad(edfadis)
 266
 267       implicit real*8 (a-h,o-z)
 268       include 'DIMENSIONS'
 269       include 'COMMON.CHAIN'
 270       include 'COMMON.DERIV'
 271       include 'COMMON.DFA'
 272
 273       double precision edfadis
 274       integer i, iatm1, iatm2,idiff
 275       double precision ckk, sckk,dist,texp
 276       double precision jix,jiy,jiz,ep,fp,scc
 277
 278       edfadis=0
 279       gdfad=0.0d0
 280
 281       do i=1, idfadis
 282
 283          iatm1=idislis(1,i)+ishiftca
 284          iatm2=idislis(2,i)+ishiftca
 285          idiff = abs(iatm1-iatm2)
 286
 287          JIX=c(1,iatm2)-c(1,iatm1)
 288          JIY=c(2,iatm2)-c(2,iatm1)
 289          JIZ=c(3,iatm2)-c(3,iatm1)
 290          DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
 291
 292          ckk=ck(idiff)
 293          sckk=sck(idiff)
 294
 295          scc = 0.0d0
 296          ep = 0.0d0
 297          fp = 0.0d0
 298
 299          do j=1,idisnum(i)
 300
 301             dd = dist-fdist(i,j)
 302             dtmp = dd*dd/ckk
 303             if (dtmp.ge.15.0d0) then
 304                texp = 0.0d0
 305             else
 306 c               texp = dfaexp( idint(dtmp*1000)+1 )/sckk
 307                 texp = exp(-dtmp)/sckk
 308             endif
 309
 310             ep=ep+sccdist(i,j)*texp
 311             fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
 312             scc=scc+sccdist(i,j)
 313 C            write(*,'(2i8,6f12.5)') i, j, dist,
 314 C     &           fdist(i,j), ep, fp, sccdist(i,j), scc
 315
 316          enddo
 317
 318          ep = -ep/scc
 319          fp = fp/scc
 320
 321
 322 c         IF(ABS(EP).lt.1.0d-20)THEN
 323 c            EP=0.0D0
 324 c         ENDIF
 325 c         IF (ABS(FP).lt.1.0d-20) THEN
 326 c            FP=0.0D0
 327 c         ENDIF
 328
 329          edfadis=edfadis+ep*dis_inc*wwdist
 330
 331          gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
 332          gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
 333          gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
 334
 335          gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
 336          gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
 337          gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
 338
 339       enddo
 340
 341       return
 342       end
 343
 344       subroutine edfat(edfator)
 345 C     DFA torsion angle
 346       implicit real*8 (a-h,o-z)
 347       include 'DIMENSIONS'
 348       include 'COMMON.CHAIN'
 349       include 'COMMON.DERIV'
 350       include 'COMMON.DFA'
 351
 352       integer i,j,ii,iii
 353       integer iatom(5)
 354       double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
 355       double precision cwidth, cwidth2
 356       PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
 357
 358       edfator= 0.0d0
 359       enephi = 0.0d0
 360       enethe = 0.0d0
 361       gdfat(:,:) = 0.0d0
 362
 363 C     START OF PHI ANGLE
 364       do i=1, idfaphi
 365
 366          aphi = 0.0d0
 367          do iii=1,5
 368           iatom(iii)=iphilis(iii,i)+ishiftca
 369          enddo
 370
 371 C     ANGLE VECTOR CALCULTION
 372          RIX=C(1,IATOM(2))-C(1,IATOM(1))
 373          RIY=C(2,IATOM(2))-C(2,IATOM(1))
 374          RIZ=C(3,IATOM(2))-C(3,IATOM(1))
 375
 376          RIPX=C(1,IATOM(3))-C(1,IATOM(2))
 377          RIPY=C(2,IATOM(3))-C(2,IATOM(2))
 378          RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
 379
 380          RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
 381          RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
 382          RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
 383
 384          RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
 385          RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
 386          RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
 387
 388          GIX=RIY*RIPZ-RIZ*RIPY
 389          GIY=RIZ*RIPX-RIX*RIPZ
 390          GIZ=RIX*RIPY-RIY*RIPX
 391
 392          GIPX=RIPY*RIPPZ-RIPZ*RIPPY
 393          GIPY=RIPZ*RIPPX-RIPX*RIPPZ
 394          GIPZ=RIPX*RIPPY-RIPY*RIPPX
 395
 396          CIPX=C(1,IATOM(3))-C(1,IATOM(1))
 397          CIPY=C(2,IATOM(3))-C(2,IATOM(1))
 398          CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
 399
 400          CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
 401          CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
 402          CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
 403
 404          CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
 405          CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
 406          CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
 407
 408          DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
 409          DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
 410          DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
 411          DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
 412
 413 C     END OF ANGLE VECTOR CALCULTION
 414
 415          TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
 416          APHI(1)=TDOT/(DGI*DRIPP)
 417          TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
 418          APHI(2)=TDOT/(DGIP*DRIP3)
 419
 420          ephi = 0.0d0
 421          tfphi1=0.0d0
 422          tfphi2=0.0d0
 423          scc=0.0d0
 424
 425          do j=1, iphinum(i)
 426             DDPS1=APHI(1)-FPHI1(i,j)
 427             DDPS2=APHI(2)-FPHI2(i,j)
 428
 429             DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
 430
 431             if (dtmp.ge.15.0d0) then
 432                ps_tmp = 0.0d0
 433             else
 434 c               ps_tmp = dfaexp(idint(dtmp*1000)+1)
 435                 ps_tmp = exp(-dtmp)
 436             endif
 437
 438             ephi=ephi+sccphi(i,j)*ps_tmp
 439
 440             tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
 441             tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
 442
 443             scc=scc+sccphi(i,j)
 444 C            write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
 445 C     &           aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
 446          ENDDO
 447
 448          ephi=-ephi/scc*phi_inc*wwangle
 449          tfphi1=tfphi1/scc*phi_inc*wwangle
 450          tfphi2=tfphi2/scc*phi_inc*wwangle
 451
 452          IF (ABS(EPHI).LT.1d-20) THEN
 453             EPHI=0.0D0
 454          ENDIF
 455          IF (ABS(TFPHI1).LT.1d-20) THEN
 456             TFPHI1=0.0D0
 457          ENDIF
 458          IF (ABS(TFPHI2).LT.1d-20) THEN
 459             TFPHI2=0.0D0
 460          ENDIF
 461
 462 C     FORCE DIRECTION CALCULATION
 463          TDX(1:5)=0.0D0
 464          TDY(1:5)=0.0D0
 465          TDZ(1:5)=0.0D0
 466
 467          DM1=1.0d0/(DGI*DRIPP)
 468
 469          GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
 470          DM2=GIRPP/(DGI**3*DRIPP)
 471          DM3=GIRPP/(DGI*DRIPP**3)
 472
 473          DM4=1.0d0/(DGIP*DRIP3)
 474
 475          GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
 476          DM5=GIRP3/(DGIP**3*DRIP3)
 477          DM6=GIRP3/(DGIP*DRIP3**3)
 478 C     FIRST ATOM BY PHI1
 479          TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
 480      &        +( GIZ* RIPY- GIY* RIPZ)*DM2
 481          TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
 482      &        +( GIX* RIPZ- GIZ* RIPX)*DM2
 483          TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
 484      &        +( GIY* RIPX- GIX* RIPY)*DM2
 485          TDX(1)=TDX(1)*TFPHI1
 486          TDY(1)=TDY(1)*TFPHI1
 487          TDZ(1)=TDZ(1)*TFPHI1
 488 C     SECOND ATOM BY PHI1
 489          TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
 490      &        -(CIPY*GIZ-CIPZ*GIY)*DM2
 491          TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
 492      &        -(CIPZ*GIX-CIPX*GIZ)*DM2
 493          TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
 494      &        -(CIPX*GIY-CIPY*GIX)*DM2
 495          TDX(2)=TDX(2)*TFPHI1
 496          TDY(2)=TDY(2)*TFPHI1
 497          TDZ(2)=TDZ(2)*TFPHI1
 498 C     SECOND ATOM BY PHI2
 499          TDX(2)=TDX(2)+
 500      &        ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
 501      &        +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
 502          TDY(2)=TDY(2)+
 503      &        ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
 504      &        +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
 505          TDZ(2)=TDZ(2)+
 506      &        ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
 507      &        +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
 508 C     THIRD ATOM BY PHI1
 509          TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
 510      &        -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
 511          TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
 512      &        -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
 513          TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
 514      &        -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
 515          TDX(3)=TDX(3)*TFPHI1
 516          TDY(3)=TDY(3)*TFPHI1
 517          TDZ(3)=TDZ(3)*TFPHI1
 518 C     THIRD ATOM BY PHI2
 519          TDX(3)=TDX(3)+
 520      &        ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
 521      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
 522          TDY(3)=TDY(3)+
 523      &        ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
 524      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
 525          TDZ(3)=TDZ(3)+
 526      &        ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
 527      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
 528 C     FOURTH ATOM BY PHI1
 529          TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
 530          TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
 531          TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
 532 C     FOURTH ATOM BY PHI2
 533          TDX(4)=TDX(4)+
 534      &        ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
 535      &        -( GIPY*RIPZ-RIPY*GIPZ)*DM5
 536      &        + RIP3X*DM6)*TFPHI2
 537          TDY(4)=TDY(4)+
 538      &        ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
 539      &        -( GIPZ*RIPX-RIPZ*GIPX)*DM5
 540      &        + RIP3Y*DM6)*TFPHI2
 541          TDZ(4)=TDZ(4)+
 542      &        ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
 543      &        -( GIPX*RIPY-RIPX*GIPY)*DM5
 544      &        + RIP3Z*DM6)*TFPHI2
 545 C     FIFTH ATOM BY PHI2
 546          TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
 547          TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
 548          TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
 549 C     END OF FORCE DIRECTION
 550 c     force calcuation
 551          DO II=1,5
 552             gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
 553             gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
 554             gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
 555          ENDDO
 556 c     energy calculation
 557          enephi = enephi + ephi
 558 c     end of single assignment statement
 559       ENDDO
 560 C     END OF PHI RESTRAINT
 561
 562 C     START OF THETA ANGLE
 563       do i=1, idfathe
 564
 565          athe = 0.0d0
 566          do iii=1,5
 567           iatom(iii)=ithelis(iii,i)+ishiftca
 568          enddo
 569
 570
 571 C     ANGLE VECTOR CALCULTION
 572          RIX=C(1,IATOM(2))-C(1,IATOM(1))
 573          RIY=C(2,IATOM(2))-C(2,IATOM(1))
 574          RIZ=C(3,IATOM(2))-C(3,IATOM(1))
 575
 576          RIPX=C(1,IATOM(3))-C(1,IATOM(2))
 577          RIPY=C(2,IATOM(3))-C(2,IATOM(2))
 578          RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
 579
 580          RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
 581          RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
 582          RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
 583
 584          RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
 585          RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
 586          RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
 587
 588          GIX=RIY*RIPZ-RIZ*RIPY
 589          GIY=RIZ*RIPX-RIX*RIPZ
 590          GIZ=RIX*RIPY-RIY*RIPX
 591
 592          GIPX=RIPY*RIPPZ-RIPZ*RIPPY
 593          GIPY=RIPZ*RIPPX-RIPX*RIPPZ
 594          GIPZ=RIPX*RIPPY-RIPY*RIPPX
 595
 596          GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
 597          GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
 598          GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
 599
 600          CIPX=C(1,IATOM(3))-C(1,IATOM(1))
 601          CIPY=C(2,IATOM(3))-C(2,IATOM(1))
 602          CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
 603
 604          CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
 605          CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
 606          CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
 607
 608          CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
 609          CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
 610          CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
 611
 612          DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
 613          DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
 614          DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
 615          DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
 616          DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
 617 C     END OF ANGLE VECTOR CALCULTION
 618
 619          TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
 620          ATHE(1)=TDOT/(DGI*DGIP)
 621          TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
 622          ATHE(2)=TDOT/(DGIP*DGIPP)
 623
 624          ETHE=0.0D0
 625          TFTHE1=0.0D0
 626          TFTHE2=0.0D0
 627          SCC=0.0D0
 628          TH_TMP=0.0d0
 629
 630          do j=1,ithenum(i)
 631             ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
 632             ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
 633             dtmp= (ddth1**2+ddth2**2)/cwidth2
 634             if ( dtmp .ge. 15.0d0) then
 635                th_tmp = 0.0d0
 636             else
 637 c               th_tmp = dfaexp ( idint(dtmp*1000)+1 )
 638                th_tmp = exp(-dtmp)
 639             end if
 640
 641             ethe=ethe+sccthe(i,j)*th_tmp
 642
 643             tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
 644             tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
 645             scc=scc+sccthe(i,j)
 646 C            write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
 647 C     &           athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
 648          enddo
 649
 650          ethe=-ethe/scc*the_inc*wwangle
 651          tfthe1=tfthe1/scc*the_inc*wwangle
 652          tfthe2=tfthe2/scc*the_inc*wwangle
 653
 654          IF (ABS(ETHE).LT.TENM20) THEN
 655             ETHE=0.0D0
 656          ENDIF
 657          IF (ABS(TFTHE1).LT.TENM20) THEN
 658             TFTHE1=0.0D0
 659          ENDIF
 660          IF (ABS(TFTHE2).LT.TENM20) THEN
 661             TFTHE2=0.0D0
 662          ENDIF
 663
 664          TDX(1:5)=0.0D0
 665          TDY(1:5)=0.0D0
 666          TDZ(1:5)=0.0D0
 667
 668          DM1=1.0d0/(DGI*DGIP)
 669          DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
 670          DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
 671
 672          DM4=1.0d0/(DGIP*DGIPP)
 673          DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
 674          DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
 675
 676 C     FIRST ATOM BY THETA1
 677          TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
 678      &        -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
 679          TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
 680      &        -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
 681          TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
 682      &        -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
 683 C     SECOND ATOM BY THETA1
 684          TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
 685      &        -(CIPY*GIZ-CIPZ*GIY)*DM2
 686      &        +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
 687          TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
 688      &        -(CIPZ*GIX-CIPX*GIZ)*DM2
 689      &        +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
 690          TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
 691      &        -(CIPX*GIY-CIPY*GIX)*DM2
 692      &        +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
 693 C     SECOND ATOM BY THETA2
 694          TDX(2)=TDX(2)+
 695      &        ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
 696      &        -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
 697          TDY(2)=TDY(2)+
 698      &        ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
 699      &        -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
 700          TDZ(2)=TDZ(2)+
 701      &        ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
 702      &        -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
 703 C     THIRD ATOM BY THETA1
 704          TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
 705      &        -(GIY*RIZ-GIZ*RIY)*DM2
 706      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
 707          TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
 708      &        -(GIZ*RIX-GIX*RIZ)*DM2
 709      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
 710          TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
 711      &        -(GIX*RIY-GIY*RIX)*DM2
 712      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
 713 C     THIRD ATOM BY THETA2
 714          TDX(3)=TDX(3)+
 715      &        ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
 716      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
 717      &        +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
 718          TDY(3)=TDY(3)+
 719      &        ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
 720      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
 721      &        +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
 722          TDZ(3)=TDZ(3)+
 723      &        ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
 724      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM5
 725      &        +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
 726 C     FOURTH ATOM BY THETA1
 727          TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
 728      &        -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
 729          TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
 730      &        -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
 731          TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
 732      &        -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
 733 C     FOURTH ATOM BY THETA2
 734          TDX(4)=TDX(4)+
 735      &        ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
 736      &        -(GIPY*RIPZ-GIPZ*RIPY)*DM5
 737      &        -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
 738          TDY(4)=TDY(4)+
 739      &        ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
 740      &        -(GIPZ*RIPX-GIPX*RIPZ)*DM5
 741      &        -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
 742          TDZ(4)=TDZ(4)+
 743      &        ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
 744      &        -(GIPX*RIPY-GIPY*RIPX)*DM5
 745      &        -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
 746 C     FIFTH ATOM BY THETA2
 747          TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
 748      &        -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
 749          TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
 750      &        -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
 751          TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
 752      &        -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
 753 C     !! END OF FORCE DIRECTION!!!!
 754          DO II=1,5
 755             gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
 756             gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
 757             gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
 758          ENDDO
 759 C     energy calculation
 760          enethe = enethe + ethe
 761       ENDDO
 762
 763       edfator = enephi + enethe
 764
 765       RETURN
 766       END
 767
 768       subroutine edfan(edfanei)
 769 C     DFA neighboring CA restraint
 770       implicit real*8 (a-h,o-z)
 771       include 'DIMENSIONS'
 772       include 'COMMON.CHAIN'
 773       include 'COMMON.DERIV'
 774       include 'COMMON.DFA'
 775
 776       integer i,j,imin
 777       integer kshnum, n1atom
 778
 779       double precision enenei,tmp_n
 780       double precision pai,hpai
 781       double precision jix,jiy,jiz,ndiff,snorm_nei
 782       double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
 783       double precision dr,dr2,half,ntmp,dtmp
 784
 785       parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
 786       parameter(pai=3.14159265358979323846D0)
 787       parameter(hpai=1.5707963267948966D0)
 788       parameter(snorm_nei=0.886226925452758D0)
 789
 790       edfanei = 0.0d0
 791       enenei  = 0.0d0
 792       gdfan   = 0.0d0
 793
 794 c      print*, 's1:', s1(:)
 795 c      print*, 's2:', s2(:)
 796
 797       do i=1, idfanei
 798
 799          kshnum=kshell(i)
 800          n1atom=ineilis(i)+ishiftca
 801 C         write(*,*) 'kshnum,n1atom:', kshnum, n1atom
 802
 803          tmp_n=0.0d0
 804          ftmp=0.0d0
 805          dnei=0.0d0
 806          dist=0.0d0
 807          t1dx=0.0d0
 808          t1dy=0.0d0
 809          t1dz=0.0d0
 810          t2dx=0.0d0
 811          t2dy=0.0d0
 812          t2dz=0.0d0
 813
 814          do j = ishiftca+1, ilastca
 815
 816             if (n1atom.eq.j) cycle
 817
 818             jix=c(1,j)-c(1,n1atom)
 819             jiy=c(2,j)-c(2,n1atom)
 820             jiz=c(3,j)-c(3,n1atom)
 821             dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
 822
 823 c            write(*,*) n1atom, j, dist
 824
 825             if(kshnum.ne.1)then
 826                if (dist.lt.s1(kshnum).and.
 827      &              dist.gt.s2(kshnum-1)) then
 828
 829                   tmp_n=tmp_n+1.0d0
 830
 831 c                  write(*,*) 'case1:',tmp_n
 832
 833 cc                  t1dx=t1dx+0.0d0
 834 cc                  t1dy=t1dy+0.0d0
 835 cc                  t1dz=t1dz+0.0d0
 836                   t2dx(j)=0.0d0
 837                   t2dy(j)=0.0d0
 838                   t2dz(j)=0.0d0
 839
 840                elseif(dist.ge.s1(kshnum).and.
 841      &                 dist.le.s2(kshnum)) then
 842
 843                   dnei=(dist-s1(kshnum))/dr2*pai
 844                   tmp_n=tmp_n + half*(1+cos(dnei))
 845 c                  write(*,*) 'case2:',tmp_n
 846                   ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
 847 c     center atom
 848                   t1dx=t1dx+jix*ftmp
 849                   t1dy=t1dy+jiy*ftmp
 850                   t1dz=t1dz+jiz*ftmp
 851 c     neighbor atoms
 852                   t2dx(j)=-jix*ftmp
 853                   t2dy(j)=-jiy*ftmp
 854                   t2dz(j)=-jiz*ftmp
 855 c
 856                elseif(dist.ge.s1(kshnum-1).and.
 857      &                 dist.le.s2(kshnum-1)) then
 858                   dnei=(dist-s1(kshnum-1))/dr2*pai
 859                   tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
 860 c                  write(*,*) 'case3:',tmp_n
 861                   ftmp = hpai*sin(dnei)/dr2/dist
 862 c     center atom
 863                   t1dx=t1dx+jix*ftmp
 864                   t1dy=t1dy+jiy*ftmp
 865                   t1dz=t1dz+jiz*ftmp
 866 c     neighbor atoms
 867                   t2dx(j)=-jix*ftmp
 868                   t2dy(j)=-jiy*ftmp
 869                   t2dz(j)=-jiz*ftmp
 870
 871                endif
 872
 873             elseif(kshnum.eq.1) then
 874
 875                if(dist.lt.s1(kshnum))then
 876
 877                   tmp_n=tmp_n+1.0d0
 878 c                  write(*,*) 'case4:',tmp_n
 879 cc                  t1dx=t1dx+0.0d0
 880 cc                  t1dy=t1dy+0.0d0
 881 cc                  t1dz=t1dz+0.0d0
 882                   t2dx(j)=0.0d0
 883                   t2dy(j)=0.0d0
 884                   t2dz(j)=0.0d0
 885
 886                elseif(dist.ge.s1(kshnum).and.
 887      &                 dist.le.s2(kshnum))then
 888
 889                   dnei=(dist-s1(kshnum))/dr2*pai
 890                   tmp_n=tmp_n + half*(1+cos(dnei))
 891 c                  write(*,*) 'case5:',tmp_n
 892                   ftmp = -hpai*sin(dnei)/dr2/dist
 893 c     center atom
 894                   t1dx=t1dx+jix*ftmp
 895                   t1dy=t1dy+jiy*ftmp
 896                   t1dz=t1dz+jiz*ftmp
 897 c     neighbor atoms
 898                   t2dx(j)=-jix*ftmp
 899                   t2dy(j)=-jiy*ftmp
 900                   t2dz(j)=-jiz*ftmp
 901
 902                endif
 903             endif
 904          enddo
 905
 906          scc=0.0d0
 907          enei=0.0d0
 908          tmp_fnei=0.0d0
 909          ndiff=0.0d0
 910
 911          do imin=1,ineinum(i)
 912
 913             ndiff = tmp_n-fnei(i,imin)
 914             dtmp  = ndiff*ndiff
 915
 916             if (dtmp.ge.15.0d0) then
 917                ntmp = 0.0d0
 918             else
 919 c               ntmp = dfaexp( idint(dtmp*1000) + 1 )
 920                 ntmp = exp(-dtmp)
 921             end if
 922
 923             enei=enei+sccnei(i,imin)*ntmp
 924             tmp_fnei=tmp_fnei-
 925      &           sccnei(i,imin)*ntmp*ndiff*2.0d0
 926             scc=scc+sccnei(i,imin)
 927
 928 c            write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
 929 c     &           fnei(i,imin),sccnei(i,imin),enei,scc
 930          enddo
 931
 932          enei=-enei/scc*snorm_nei*nei_inc*wwnei
 933          tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
 934
 935 c         if (abs(enei).lt.1.0d-20)then
 936 c            enei=0.0d0
 937 c         endif
 938 c         if (abs(tmp_fnei).lt.1.0d-20) then
 939 c            tmp_fnei=0.0d0
 940 c         endif
 941
 942 c     force calculation
 943          t1dx=t1dx*tmp_fnei
 944          t1dy=t1dy*tmp_fnei
 945          t1dz=t1dz*tmp_fnei
 946
 947          do j=ishiftca+1,ilastca
 948             t2dx(j)=t2dx(j)*tmp_fnei
 949             t2dy(j)=t2dy(j)*tmp_fnei
 950             t2dz(j)=t2dz(j)*tmp_fnei
 951          enddo
 952
 953          gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
 954          gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
 955          gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
 956
 957          do j=ishiftca+1,ilastca
 958             gdfan(1,j)=gdfan(1,j)+t2dx(j)
 959             gdfan(2,j)=gdfan(2,j)+t2dy(j)
 960             gdfan(3,j)=gdfan(3,j)+t2dz(j)
 961          enddo
 962 c     energy calculation
 963
 964          enenei=enenei+enei
 965
 966       enddo
 967
 968       edfanei=enenei
 969
 970       return
 971       end
 972
 973       subroutine edfab(edfabeta)
 974
 975       implicit real*8 (a-h,o-z)
 976
 977       include 'DIMENSIONS'
 978       include 'COMMON.CHAIN'
 979       include 'COMMON.DERIV'
 980       include 'COMMON.DFA'
 981
 982       real*8 PAI
 983       parameter(PAI=3.14159265358979323846D0)
 984       parameter (maxca=800)
 985 C     sheet variables
 986       real*8 bx(maxres),by(maxres),bz(maxres)
 987       real*8 vbet(maxres,maxres)
 988       real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
 989       real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
 990       real*8 vbeta,vbetp,vbetm
 991       real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
 992      &     c00,s00,ulnex,dnex
 993       real*8 dp45,dm45,w_beta
 994
 995       real*8 cph(maxca),cth(maxca)
 996       real*8 atx(maxca),aty(maxca),atz(maxca)
 997       real*8 atmx(maxca),atmy(maxca),atmz(maxca)
 998       real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
 999       real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1000       real*8 sth(maxca)
1001       real*8 astx(maxca),asty(maxca),astz(maxca)
1002       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1003       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1004       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1005
1006       real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
1007      & astxnum(maxca),astynum(maxca),astznum(maxca),
1008      & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
1009      & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
1010      & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
1011      & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
1012      & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
1013      & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
1014      & cth_orig(maxca),sth_orig(maxca)
1015
1016       common /sheca/     bx,by,bz
1017       common /shee/      vbeta,vbet,vbetp,vbetm
1018       common /shetf/     shetfx,shetfy,shetfz
1019       common /shef/      shefx, shefy, shefz
1020       common /sheparm/   dca,dlhb,ulhb,dshe,dldhb,uldhb,
1021      &                   c00,s00,ulnex,dnex
1022       common /sheconst/  dp45,dm45,w_beta
1023
1024       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1025      $     atmmz,atm3x,atm3y,atm3z
1026       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1027      $     astmmz,astm3x,astm3y,astm3z
1028
1029       common /coscos/   cph,cth
1030       common /sinsin/ sth
1031
1032 C     End of sheet variables
1033
1034       integer i,j
1035       double precision enebet
1036
1037       enebet=0.0d0
1038 c      bx=0.0d0;by=0.0d0;bz=0.0d0
1039 c      shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1040
1041       gdfab=0.0d0
1042
1043       do i=ishiftca+1,ilastca
1044          bx(i-ishiftca)=c(1,i)
1045          by(i-ishiftca)=c(2,i)
1046          bz(i-ishiftca)=c(3,i)
1047       enddo
1048
1049 c      do i=1,ilastca-ishiftca
1050 c         read(99,*) bx(i),by(i),bz(i)
1051 c      enddo
1052 c      close(99)
1053
1054       dca=0.25d0**2
1055       dshe=0.3d0**2
1056       ULHB=5.0D0
1057       ULDHB=5.0D0
1058       ULNEX=COS(60.0D0/180.0D0*PAI)
1059
1060       DLHB=1.0D0
1061       DLDHB=1.0D0
1062
1063       DNEX=0.3D0**2
1064
1065       C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1066       S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1067
1068       W_BETA=0.5D0
1069       DP45=W_BETA
1070       DM45=W_BETA
1071
1072 C     END OF INITIALIZATION
1073
1074       nca=ilastca-ishiftca
1075
1076       call angvectors(nca)
1077       call sheetforce(nca,wshet)
1078
1079 c     end of sheet energy and force
1080
1081       do j=1,nca
1082          shetfx(j)=shetfx(j)*beta_inc
1083          shetfy(j)=shetfy(j)*beta_inc
1084          shetfz(j)=shetfz(j)*beta_inc
1085 c         write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1086       enddo
1087
1088       vbeta=vbeta*beta_inc
1089       enebet=vbeta
1090       edfabeta=enebet
1091
1092       do j=1,nca
1093          gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1094          gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1095          gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1096       enddo
1097
1098 #ifdef DEBUG1
1099       do j=1,nca
1100         write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1101       enddo
1102
1103
1104       gdfab=0
1105       dinc=0.001
1106       do j=1,nca
1107         cth_orig(j)=cth(j)
1108         sth_orig(j)=sth(j)
1109       enddo
1110
1111       do j=1,nca
1112
1113        bx(j)=bx(j)+dinc
1114        call angvectors(nca)
1115        bx(j)=bx(j)-2*dinc
1116        call angvectors(nca)
1117        atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1118        astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1119        if (j.gt.1) then
1120        atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1121        astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1122        endif
1123        if (j.gt.2) then
1124        atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1125        astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1126        endif
1127        if (j.gt.3) then
1128        atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1129        astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1130        endif
1131        bx(j)=bx(j)+dinc
1132        by(j)=by(j)+dinc
1133        call angvectors(nca)
1134        by(j)=by(j)-2*dinc
1135        call angvectors(nca)
1136        by(j)=by(j)+dinc
1137        atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1138        astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1139        if (j.gt.1) then
1140        atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1141        astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1142        endif
1143        if (j.gt.2) then
1144        atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1145        astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1146        endif
1147        if (j.gt.3) then
1148        atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1149        astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1150        endif
1151
1152        bz(j)=bz(j)+dinc
1153        call angvectors(nca)
1154        bz(j)=bz(j)-2*dinc
1155        call angvectors(nca)
1156        bz(j)=bz(j)+dinc
1157
1158        atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1159        astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1160        if (j.gt.1) then
1161        atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1162        astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1163        endif
1164        if (j.gt.2) then
1165        atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1166        astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1167        endif
1168        if (j.gt.3) then
1169        atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1170        astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1171        endif
1172
1173       enddo
1174
1175       do i=1,nca
1176         write (*,'(2i5,a2,6f10.5)')
1177      &  i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1178      &          astxnum(i),astx(i),astxnum(i)/astx(i),
1179      &  i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1180      &          astynum(i),asty(i),astynum(i)/asty(i),
1181      &  i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1182      &          astznum(i),astz(i),astznum(i)/astz(i),
1183      &  i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1184      &          astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1185      &  i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1186      &          astmynum(i),astmy(i),astmynum(i)/astmy(i),
1187      &  i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1188      &          astmznum(i),astmz(i),astmznum(i)/astmz(i),
1189      &  i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1190      &          astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1191      &  i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1192      &          astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1193      &  i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1194      &          astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1195      &  i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1196      &          astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1197      &  i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1198      &          astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1199      &  i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1200      &          astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1201      &  i,0," ",cth_orig(i),sth_orig(i)
1202       enddo
1203
1204
1205       gdfab=0
1206       dinc=0.001
1207
1208       do j=1,nca
1209
1210        bx(j)=bx(j)+dinc
1211        call angvectors(nca)
1212        call sheetforce(nca,wshet)
1213        vbeta1=vbeta*beta_inc
1214        bx(j)=bx(j)-2*dinc
1215        call angvectors(nca)
1216        call sheetforce(nca,wshet)
1217        vbeta2=vbeta*beta_inc
1218        gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1219        bx(j)=bx(j)+dinc
1220
1221        by(j)=by(j)+dinc
1222        call angvectors(nca)
1223        call sheetforce(nca,wshet)
1224        vbeta1=vbeta*beta_inc
1225        by(j)=by(j)-2*dinc
1226        call angvectors(nca)
1227        call sheetforce(nca,wshet)
1228        vbeta2=vbeta*beta_inc
1229        gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1230        by(j)=by(j)+dinc
1231
1232        bz(j)=bz(j)+dinc
1233        call angvectors(nca)
1234        call sheetforce(nca,wshet)
1235        vbeta1=vbeta*beta_inc
1236        bz(j)=bz(j)-2*dinc
1237        call angvectors(nca)
1238        call sheetforce(nca,wshet)
1239        vbeta2=vbeta*beta_inc
1240        gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1241        bz(j)=bz(j)+dinc
1242
1243
1244       enddo
1245
1246
1247       call angvectors(nca)
1248       call sheetforce(nca,wshet)
1249       do j=1,nca
1250          shetfx(j)=shetfx(j)*beta_inc
1251          shetfy(j)=shetfy(j)*beta_inc
1252          shetfz(j)=shetfz(j)*beta_inc
1253       enddo
1254
1255
1256       write(*,*) 'xyz analytical and numerical gradient'
1257       do j=1,nca
1258         write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1259      &                   ,(-gdfab(i,j),i=1,3)
1260       enddo
1261
1262       do j=1,nca
1263         write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1264      &                                  shetfy(j)/gdfab(2,j),
1265      &                                  shetfz(j)/gdfab(3,j)
1266       enddo
1267
1268       stop
1269 #endif
1270
1271       return
1272       end
1273 C-------------------------------------------------------------------------------
1274       subroutine angvectors(nca)
1275 c      implicit real*4(a-h,o-z)
1276       implicit none
1277       integer nca
1278       integer maxca
1279       parameter(maxca=800)
1280       real*8   pai,zero
1281       parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1282
1283       real*8   bx(maxca),by(maxca),bz(maxca)
1284       real*8   dis(maxca,maxca)
1285       real*8   apx(maxca),apy(maxca),apz(maxca)
1286       real*8   apmx(maxca),apmy(maxca),apmz(maxca)
1287       real*8   apmmx(maxca),apmmy(maxca),apmmz(maxca)
1288       real*8   apm3x(maxca),apm3y(maxca),apm3z(maxca)
1289       real*8   atx(maxca),aty(maxca),atz(maxca)
1290       real*8   atmx(maxca),atmy(maxca),atmz(maxca)
1291       real*8   atmmx(maxca),atmmy(maxca),atmmz(maxca)
1292       real*8   atm3x(maxca),atm3y(maxca),atm3z(maxca)
1293       real*8   astx(maxca),asty(maxca),astz(maxca)
1294       real*8   astmx(maxca),astmy(maxca),astmz(maxca)
1295       real*8   astmmx(maxca),astmmy(maxca),astmmz(maxca)
1296       real*8   astm3x(maxca),astm3y(maxca),astm3z(maxca)
1297       real*8   sth(maxca)
1298       real*8   cph(maxca),cth(maxca)
1299       real*8   ulcos(maxca)
1300       real*8   p,c
1301       integer  i, ip, ipp, ip3, j
1302       real*8   rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1303       real*8   rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1304       real*8   gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1305       real*8   cix, ciy, ciz, cipx, cipy, cipz
1306       real*8   gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1307       real*8   d10, d11, d12, d13, d20, d21, d22, d23, d24
1308       real*8   d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1309       real*8   d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1310       real*8   dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1311       real*8   dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1312       real*8   g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1313       real*8   gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1314       real*8   gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1315       real*8   gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1316       real*8   gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1317       real*8   grpp,gx,gy,gz
1318       real*8   rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1319       real*8   sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1320       integer inb,nmax,iselect
1321
1322       common /sheca/   bx,by,bz
1323       common /difvec/  rx, ry, rz
1324       common /ulang/    ulcos
1325       common /phys1/   inb,nmax,iselect
1326       common /phys4/   p,c
1327       common /kyori2/  dis
1328       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1329      &     apmmz,apm3x,apm3y,apm3z
1330       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1331      &     atmmz,atm3x,atm3y,atm3z
1332       common /coscos/   cph,cth
1333       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1334      &     astmmz,astm3x,astm3y,astm3z
1335       common /sinsin/   sth
1336 C-------------------------------------------------------------------------------
1337 c      write(*,*) 'inside angvectors'
1338 C     initialize
1339       p=0.1d0
1340       c=1.0d0
1341       inb=nca
1342       cph=zero; cth=zero; sth=zero
1343       apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1344       apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1345       atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1346       atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1347       astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1348       astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1349       astm3z=zero
1350 C     end of initialize
1351 C     r[x,y,z] calc and distance calculation
1352       rx=zero;ry=zero;rz=zero
1353
1354       do i=1,inb
1355          do j=1,inb
1356             rx(i,j)=bx(j)-bx(i)
1357             ry(i,j)=by(j)-by(i)
1358             rz(i,j)=bz(j)-bz(i)
1359             dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1360 c            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1361 c            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1362 c            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1363 c            write(*,*) 'dis(i,j):',i,j,dis(i,j)
1364          enddo
1365       enddo
1366 c     end of r[x,y,z] calc
1367 C     cos calc
1368       do i=1,inb-2
1369          ip=i+1
1370          ipp=i+2
1371
1372          if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1373             ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1374      $           +rz(i,ip)*rz(ip,ipp)
1375             ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1376          endif
1377       enddo
1378 c     end of virtual bond angle
1379 c      write(*,*) 'inside angvectors1'
1380 crc       do i=1,inb-3
1381       do i=1,inb
1382          ip=i+1
1383          ipp=i+2
1384          ip3=i+3
1385          rix=bx(ip)-bx(i)
1386          riy=by(ip)-by(i)
1387          riz=bz(ip)-bz(i)
1388          ripx=bx(ipp)-bx(ip)
1389          ripy=by(ipp)-by(ip)
1390          ripz=bz(ipp)-bz(ip)
1391          rippx=bx(ip3)-bx(ipp)
1392          rippy=by(ip3)-by(ipp)
1393          rippz=bz(ip3)-bz(ipp)
1394
1395          gx=riy*ripz-riz*ripy
1396          gy=riz*ripx-rix*ripz
1397          gz=rix*ripy-riy*ripx
1398          gpx=ripy*rippz-ripz*rippy
1399          gpy=ripz*rippx-ripx*rippz
1400          gpz=ripx*rippy-ripy*rippx
1401          gpcrp_x=gpy*ripz-gpz*ripy
1402          gpcrp_y=gpz*ripx-gpx*ripz
1403          gpcrp_z=gpx*ripy-gpy*ripx
1404          d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1405          gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1406      &        -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1407
1408          if(i.ge.2) then
1409             rimx=bx(i)-bx(i-1)
1410             rimy=by(i)-by(i-1)
1411             rimz=bz(i)-bz(i-1)
1412             gmx=rimy*riz-rimz*riy
1413             gmy=rimz*rix-rimx*riz
1414             gmz=rimx*riy-rimy*rix
1415             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1416             dgm3=dgm**3
1417             ggm=gmx*gx+gmy*gy+gmz*gz
1418             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1419             drim=dis(i-1,i)
1420             drim3=drim**3
1421             gcr_x=gy*riz-gz*riy
1422             gcr_y=gz*rix-gx*riz
1423             gcr_z=gx*riy-gy*rix
1424             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1425             d_gcr3=d_gcr**3
1426             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1427      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1428          endif
1429 c         write(*,*) 'inside angvectors2'
1430          if(i.ge.3) then
1431             rimmx=bx(i-1)-bx(i-2)
1432             rimmy=by(i-1)-by(i-2)
1433             rimmz=bz(i-1)-bz(i-2)
1434             drimm=dis(i-2,i-1)
1435             gmmx=rimmy*rimz-rimmz*rimy
1436             gmmy=rimmz*rimx-rimmx*rimz
1437             gmmz=rimmx*rimy-rimmy*rimx
1438             dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1439             dgmm3=dgmm**3
1440             gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1441             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1442             gmcrim_x=gmy*rimz-gmz*rimy
1443             gmcrim_y=gmz*rimx-gmx*rimz
1444             gmcrim_z=gmx*rimy-gmy*rimx
1445             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1446             d_gmcrim3=d_gmcrim**3
1447             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1448      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1449          endif
1450
1451          if(i.ge.4) then
1452             rim3x=bx(i-2)-bx(i-3)
1453             rim3y=by(i-2)-by(i-3)
1454             rim3z=bz(i-2)-bz(i-3)
1455             g3x=rim3y*rimmz-rim3z*rimmy
1456             g3y=rim3z*rimmx-rim3x*rimmz
1457             g3z=rim3x*rimmy-rim3y*rimmx
1458             dg30=sqrt(g3x**2+g3y**2+g3z**2)
1459             g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1460             g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1461 cc**********************************************************************
1462             gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1463             gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1464             gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1465             d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1466             d_gmmcrimm3=d_gmmcrimm**3
1467             gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1468      &           -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1469          endif
1470
1471          dri=dis(i,i+1)
1472          drip=dis(i+1,i+2)
1473          dripp=dis(i+2,i+3)
1474          dri3=dri**3
1475          dg=sqrt(gx**2+gy**2+gz**2)
1476          dgp=sqrt(gpx**2+gpy**2+gpz**2)
1477          dg3=dg**3
1478
1479          ggp=gx*gpx+gy*gpy+gz*gpz
1480          grpp=gx*rippx+gy*rippy+gz*rippz
1481
1482          if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1483      &        .and.d_gpcrp.gt.0.0D0) then
1484             cph(i)=grpp/dg/dripp
1485             cth(i)=ggp/dg/dgp
1486             sth(i)=gpcrp__g/d_gpcrp/dg
1487          else
1488 c
1489             cph(i)=1.0D0
1490             cth(i)=1.0D0
1491             sth(i)=0.0D0
1492          endif
1493
1494 c         write(*,*) 'inside angvectors3'
1495
1496          if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1497      &        .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1498             d10=1.0D0/(dg*dgp)
1499             d11=ggp/(dg3*dgp)
1500             d12=1.0D0/(dg*dripp)
1501             d13=grpp/(dg3*dripp)
1502             sd10=1.0D0/(d_gpcrp*dg)
1503             sd11=gpcrp__g/(d_gpcrp*dg3)
1504          else
1505             d10=0.0D0
1506             d11=0.0D0
1507             d12=0.0D0
1508             d13=0.0D0
1509             sd10=0.0D0
1510             sd11=0.0D0
1511          endif
1512
1513          atx(i)=(ripz*gpy-ripy*gpz)*d10
1514      &        -(gy*ripz-gz*ripy)*d11
1515          aty(i)=(ripx*gpz-ripz*gpx)*d10
1516      &        -(gz*ripx-gx*ripz)*d11
1517          atz(i)=(ripy*gpx-ripx*gpy)*d10
1518      &        -(gx*ripy-gy*ripx)*d11
1519          astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1520      &        +ripy*gpy*ripx-gpx*ripz**2)
1521      &        -sd11*(gy*ripz-gz*ripy)
1522          asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1523      &        -gpy*ripx**2+gpz*ripy*ripz)
1524      &        -sd11*(-gx*ripz+gz*ripx)
1525          astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1526      &        -gpz*ripy**2+ripz*gpx*ripx)
1527      &        -sd11*(gx*ripy-gy*ripx)
1528          apx(i)=(ripz*rippy-ripy*rippz)*d12
1529      &        -(gy*ripz-gz*ripy)*d13
1530          apy(i)=(ripx*rippz-ripz*rippx)*d12
1531      &        -(gz*ripx-gx*ripz)*d13
1532          apz(i)=(ripy*rippx-ripx*rippy)*d12
1533      &        -(gx*ripy-gy*ripx)*d13
1534
1535          if(i.ge.2) then
1536             cix=bx(ip)-bx(i-1)
1537             ciy=by(ip)-by(i-1)
1538             ciz=bz(ip)-bz(i-1)
1539             cipx=bx(ipp)-bx(i)
1540             cipy=by(ipp)-by(i)
1541             cipz=bz(ipp)-bz(i)
1542             ripx=bx(ipp)-bx(ip)
1543             ripy=by(ipp)-by(ip)
1544             ripz=bz(ipp)-bz(ip)
1545             if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1546      &           .and.d_gcr3.gt.0.0D0) then
1547                d20=1.0D0/(dg*dgm)
1548                d21=ggm/(dgm3*dg)
1549                d22=ggm/(dgm*dg3)
1550                d23=1.0D0/(dgm*drip)
1551                d24=gmrp/(dgm3*drip)
1552                sd20=1.0D0/(d_gcr*dgm)
1553                sd21=gcr__gm/(d_gcr3*dgm)
1554                sd22=gcr__gm/(d_gcr*dgm3)
1555             else
1556                d20=0.0D0
1557                d21=0.0D0
1558                d22=0.0D0
1559                d23=0.0D0
1560                d24=0.0D0
1561                sd20=0.0D0
1562                sd21=0.0D0
1563                sd22=0.0D0
1564             endif
1565             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1566      &           -(ciy*gmz-ciz*gmy)*d21
1567      &           +(ripy*gz-ripz*gy)*d22
1568             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1569      &           -(ciz*gmx-cix*gmz)*d21
1570      &           +(ripz*gx-ripx*gz)*d22
1571             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1572      &           -(cix*gmy-ciy*gmx)*d21
1573      &           +(ripx*gy-ripy*gx)*d22
1574 cc**********************************************************************
1575             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1576      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1577      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1578      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1579      &           +gcr_z*(-ripz*rix+gy))
1580      &           -sd22*(-gmy*ciz+gmz*ciy)
1581
1582             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1583      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1584      &           +riz*ripz*gmy)
1585      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1586      &           -gcr_z*(ripz*riy+gx))
1587      &           -sd22*(gmx*ciz-gmz*cix)
1588
1589             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1590      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1591      &           -riz*gx*cix)
1592      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1593      &           +gcr_z*(ripy*riy+ripx*rix))
1594      &           -sd22*(-gmx*ciy+gmy*cix)
1595 cc**********************************************************************
1596             apmx(i)=(ciy*ripz-ripy*ciz)*d23
1597      &           -(ciy*gmz-ciz*gmy)*d24
1598             apmy(i)=(ciz*ripx-ripz*cix)*d23
1599      &           -(ciz*gmx-cix*gmz)*d24
1600             apmz(i)=(cix*ripy-ripx*ciy)*d23
1601      &           -(cix*gmy-ciy*gmx)*d24
1602          endif
1603
1604          if(i.ge.3) then
1605             if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1606      &           .and.d_gmcrim3.gt.0.0D0) then
1607                d30=1.0D0/(dgm*dgmm)
1608                d31=gmmgm/(dgm3*dgmm)
1609                d32=gmmgm/(dgm*dgmm3)
1610                d33=1.0D0/(dgmm*dri)
1611                d34=gmmr/(dgmm3*dri)
1612                d35=gmmr/(dgmm*dri3)
1613                sd30=1.0D0/(d_gmcrim*dgmm)
1614                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1615                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1616             else
1617                d30=0.0D0
1618                d31=0.0D0
1619                d32=0.0D0
1620                d33=0.0D0
1621                d34=0.0D0
1622                d35=0.0D0
1623                sd30=0.0D0
1624                sd31=0.0D0
1625                sd32=0.0D0
1626             endif
1627
1628 c            write(*,*) 'inside angvectors4'
1629
1630 cc**********************************************************************
1631             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1632      &           -(ciy*gmz-ciz*gmy)*d31
1633      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1634             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1635      &           -(ciz*gmx-cix*gmz)*d31
1636      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1637             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1638      &           -(cix*gmy-ciy*gmx)*d31
1639      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1640 cc**********************************************************************
1641             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1642      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1643      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1644      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1645      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1646      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1647      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1648
1649             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1650      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1651      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1652      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1653      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1654      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1655      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1656
1657             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1658      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1659      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1660      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1661      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1662      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1663      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1664 c**********************************************************************
1665             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1666      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1667      &           +rix*d35
1668             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1669      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1670      &           +riy*d35
1671             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1672      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1673      &           +riz*d35
1674          endif
1675
1676          if(i.ge.4) then
1677             if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1678      &           .and.drim3.gt.0.0D0
1679      &           .and.d_gmmcrimm3.gt.0.0D0) then
1680                d40=1.0D0/(dg30*dgmm)
1681                d41=g3gmm/(dg30*dgmm3)
1682                d42=1.0D0/(dg30*drim)
1683                d43=g3rim_/(dg30*drim3)
1684                sd40=1.0D0/(dg30*d_gmmcrimm)
1685                sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1686             else
1687                d40=0.0D0
1688                d41=0.0D0
1689                d42=0.0D0
1690                d43=0.0D0
1691                sd40=0.0D0
1692                sd41=0.0D0
1693             endif
1694             atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1695      &           -(gmmy*rimmz-gmmz*rimmy)*d41
1696             atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1697      &           -(gmmz*rimmx-gmmx*rimmz)*d41
1698             atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1699      &           -(gmmx*rimmy-gmmy*rimmx)*d41
1700 cc**********************************************************************
1701             astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1702      &           -g3z*rimmz*rimmx+rimmy**2*g3x)
1703      &           -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1704      &           -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1705
1706             astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1707      &           -rimmx*rimmy*g3x+rimmz**2*g3y)
1708      &           -sd41*(-gmmcrimm_x*rimmx*rimmy
1709      &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1710
1711 c     &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1712
1713             astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1714      &           +g3z*rimmx**2-rimmz*rimmy*g3y)
1715      &           -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1716      &           +gmmcrimm_z*(rimmy**2+rimmx**2))
1717 c**********************************************************************
1718             apm3x(i)=g3x*d42-rimx*d43
1719             apm3y(i)=g3y*d42-rimy*d43
1720             apm3z(i)=g3z*d42-rimz*d43
1721          endif
1722       enddo
1723 c*******************************************************************************
1724
1725 c      write(*,*) 'inside angvectors5'
1726
1727 c       do i=inb-2,inb
1728        do i=1,0
1729          rimx=bx(i)-bx(i-1)
1730          rimy=by(i)-by(i-1)
1731          rimz=bz(i)-bz(i-1)
1732          rimmx=bx(i-1)-bx(i-2)
1733          rimmy=by(i-1)-by(i-2)
1734          rimmz=bz(i-1)-bz(i-2)
1735          rim3x=bx(i-2)-bx(i-3)
1736          rim3y=by(i-2)-by(i-3)
1737          rim3z=bz(i-2)-bz(i-3)
1738          gmmx=rimmy*rimz-rimmz*rimy
1739          gmmy=rimmz*rimx-rimmx*rimz
1740          gmmz=rimmx*rimy-rimmy*rimx
1741          g3x=rim3y*rimmz-rim3z*rimmy
1742          g3y=rim3z*rimmx-rim3x*rimmz
1743          g3z=rim3x*rimmy-rim3y*rimmx
1744
1745          dg30=sqrt(g3x**2+g3y**2+g3z**2)
1746          g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1747          dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1748          dgmm3=dgmm**3
1749          drim=dis(i-1,i)
1750          drimm=dis(i-2,i-1)
1751          drim3=drim**3
1752          g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1753 cc**********************************************************************
1754          gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1755          gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1756          gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1757          d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1758          d_gmmcrimm3=d_gmmcrimm**3
1759          gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1760      &        -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1761
1762          if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1763      &        .and.drim3.gt.0.0D0
1764      &        .and.d_gmmcrimm3.gt.0.0D0) then
1765             d40=1.0D0/(dg30*dgmm)
1766             d41=g3gmm/(dg30*dgmm3)
1767             d42=1.0D0/(dg30*drim)
1768             d43=g3rim_/(dg30*drim3)
1769             sd40=1.0D0/(dg30*d_gmmcrimm)
1770             sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1771          else
1772             d40=0.0D0
1773             d41=0.0D0
1774             d42=0.0D0
1775             d43=0.0D0
1776             sd40=0.0D0
1777             sd41=0.0D0
1778          endif
1779          atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1780      &        -(gmmy*rimmz-gmmz*rimmy)*d41
1781          atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1782      &        -(gmmz*rimmx-gmmx*rimmz)*d41
1783          atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1784      &        -(gmmx*rimmy-gmmy*rimmx)*d41
1785 cc**********************************************************************
1786          astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1787      &        -g3z*rimmz*rimmx+rimmy**2*g3x)
1788      &        -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1789      &        -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1790
1791          astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1792      &        -rimmx*rimmy*g3x+rimmz**2*g3y)
1793      &        -sd41*(-gmmcrimm_x*rimmx*rimmy
1794      &        +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1795
1796          astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1797      &        +g3z*rimmx**2-rimmz*rimmy*g3y)
1798      &        -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1799      &        +gmmcrimm_z*(rimmy**2+rimmx**2))
1800 cc**********************************************************************
1801          apm3x(i)=g3x*d42-rimx*d43
1802          apm3y(i)=g3y*d42-rimy*d43
1803          apm3z(i)=g3z*d42-rimz*d43
1804
1805          if(i.le.inb-1) then
1806             ip=i+1
1807             rix=bx(ip)-bx(i)
1808             riy=by(ip)-by(i)
1809             riz=bz(ip)-bz(i)
1810             cix=bx(ip)-bx(i-1)
1811             ciy=by(ip)-by(i-1)
1812             ciz=bz(ip)-bz(i-1)
1813             gmx=rimy*riz-rimz*riy
1814             gmy=rimz*rix-rimx*riz
1815             gmz=rimx*riy-rimy*rix
1816             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1817             dgm3=dgm**3
1818             dri=dis(i,i+1)
1819             dri3=dri**3
1820             gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1821             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1822             gmcrim_x=gmy*rimz-gmz*rimy
1823             gmcrim_y=gmz*rimx-gmx*rimz
1824             gmcrim_z=gmx*rimy-gmy*rimx
1825             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1826             d_gmcrim3=d_gmcrim**3
1827             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1828      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1829
1830             if(dgm3.gt.0.0D0.and.
1831      &           dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1832      &           .and.d_gmcrim3.gt.0.0D0) then
1833                d30=1.0D0/(dgm*dgmm)
1834                d31=gmmgm/(dgm3*dgmm)
1835                d32=gmmgm/(dgm*dgmm3)
1836                d33=1.0D0/(dgmm*dri)
1837                d34=gmmr/(dgmm3*dri)
1838                d35=gmmr/(dgmm*dri3)
1839                sd30=1.0D0/(d_gmcrim*dgmm)
1840                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1841                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1842
1843             else
1844                d30=0.0D0
1845                d31=0.0D0
1846                d32=0.0D0
1847                d33=0.0D0
1848                d34=0.0D0
1849                d35=0.0D0
1850                sd30=0.0D0
1851                sd31=0.0D0
1852                sd32=0.0D0
1853             endif
1854 cc**********************************************************************
1855             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1856      &           -(ciy*gmz-ciz*gmy)*d31
1857      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1858             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1859      &           -(ciz*gmx-cix*gmz)*d31
1860      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1861             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1862      &           -(cix*gmy-ciy*gmx)*d31
1863      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1864 cc**********************************************************************
1865             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1866      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1867      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1868      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1869      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1870      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1871      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1872
1873             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1874      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1875      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1876      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1877      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1878      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1879      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1880
1881             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1882      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1883      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1884      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1885      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1886      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1887      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1888 cc**********************************************************************
1889             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1890      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1891      &           +rix*d35
1892             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1893      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1894      &           +riy*d35
1895             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1896      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1897      &           +riz*d35
1898          endif
1899
1900 c         write(*,*) 'inside angvectors6'
1901
1902          if(i.eq.inb-2) then
1903             ipp=i+2
1904             ripx=bx(ipp)-bx(ip)
1905             ripy=by(ipp)-by(ip)
1906             ripz=bz(ipp)-bz(ip)
1907             cipx=bx(ipp)-bx(i)
1908             cipy=by(ipp)-by(i)
1909             cipz=bz(ipp)-bz(i)
1910             gx=riy*ripz-riz*ripy
1911             gy=riz*ripx-rix*ripz
1912             gz=rix*ripy-riy*ripx
1913             ggm=gmx*gx+gmy*gy+gmz*gz
1914             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1915             dg=sqrt(gx**2+gy**2+gz**2)
1916             dg3=dg**3
1917             drip=dis(i+1,i+2)
1918             gcr_x=gy*riz-gz*riy
1919             gcr_y=gz*rix-gx*riz
1920             gcr_z=gx*riy-gy*rix
1921             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1922             d_gcr3=d_gcr**3
1923             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1924      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1925             if(dgm3.gt.0.0D0.and.
1926      &           dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1927      &           ) then
1928                d20=1.0D0/(dg*dgm)
1929                d21=ggm/(dgm3*dg)
1930                d22=ggm/(dgm*dg3)
1931                d23=1.0D0/(dgm*drip)
1932                d24=gmrp/(dgm3*drip)
1933                sd20=1.0D0/(d_gcr*dgm)
1934                sd21=gcr__gm/(d_gcr3*dgm)
1935                sd22=gcr__gm/(d_gcr*dgm3)
1936             else
1937                d20=0.0D0
1938                d21=0.0D0
1939                d22=0.0D0
1940                d23=0.0D0
1941                d24=0.0D0
1942                sd20=0.0D0
1943                sd21=0.0D0
1944                sd22=0.0D0
1945             endif
1946             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1947      &           -(ciy*gmz-ciz*gmy)*d21
1948      &           +(ripy*gz-ripz*gy)*d22
1949             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1950      &           -(ciz*gmx-cix*gmz)*d21
1951      &           +(ripz*gx-ripx*gz)*d22
1952             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1953      &           -(cix*gmy-ciy*gmx)*d21
1954      &           +(ripx*gy-ripy*gx)*d22
1955 cc**********************************************************************
1956             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1957      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1958      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1959      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1960      &           +gcr_z*(-ripz*rix+gy))
1961      &           -sd22*(-gmy*ciz+gmz*ciy)
1962
1963             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1964      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1965      &           +riz*ripz*gmy)
1966      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1967      &           -gcr_z*(ripz*riy+gx))
1968      &           -sd22*(gmx*ciz-gmz*cix)
1969
1970             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1971      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1972      &           -riz*gx*cix)
1973      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1974      &           +gcr_z*(ripy*riy+ripx*rix))
1975      &           -sd22*(-gmx*ciy+gmy*cix)
1976 cc**********************************************************************
1977 c
1978             apmx(i)=(ciy*ripz-ripy*ciz)*d23
1979      &           -(ciy*gmz-ciz*gmy)*d24
1980             apmy(i)=(ciz*ripx-ripz*cix)*d23
1981      &           -(ciz*gmx-cix*gmz)*d24
1982             apmz(i)=(cix*ripy-ripx*ciy)*d23
1983      &           -(cix*gmy-ciy*gmx)*d24
1984
1985          endif
1986       enddo
1987
1988       return
1989       end
1990 c     END of angvectors
1991 c-------------------------------------------------------------------------------
1992 C---------------------------------------------------------------------------------
1993       subroutine sheetforce(nca,wshet)
1994       implicit none
1995 C     JYLEE
1996 c     this should be matched with dfa.fcm
1997       integer maxca
1998       parameter(maxca=800)
1999 cc**********************************************************************
2000       integer nca
2001       integer i,k
2002       integer inb,nmax,iselect
2003
2004 c      real*8 dfaexp(15001)
2005
2006       real*8 vbeta,vbetp,vbetm
2007       real*8 shefx(maxca,12)
2008       real*8 shefy(maxca,12),shefz(maxca,12)
2009       real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2010       real*8 vbet(maxca,maxca)
2011       real*8 wshet(maxca,maxca)
2012       real*8 bx(maxca),by(maxca),bz(maxca)
2013
2014       common /sheca/  bx,by,bz
2015       common /phys1/  inb,nmax,iselect
2016       common /shef/   shefx,shefy,shefz
2017       common /shee/   vbeta,vbet,vbetp,vbetm
2018       common /shetf/  shetfx,shetfy,shetfz
2019
2020       inb=nca
2021       do i=1,inb
2022          shetfx(i)=0.0D0
2023          shetfy(i)=0.0D0
2024          shetfz(i)=0.0D0
2025       enddo
2026
2027       do k=1,12
2028          do i=1,inb
2029             shefx(i,k)=0.0D0
2030             shefy(i,k)=0.0D0
2031             shefz(i,k)=0.0D0
2032          enddo
2033       enddo
2034
2035       call sheetene(nca,wshet)
2036       call sheetforce1
2037
2038  887  format(a,1x,i6,3x,f12.8)
2039  888  format(a,1x,i4,1x,i4,3x,f12.8)
2040  889  format(a,1x,i4,3x,f12.8)
2041       !write(2,*) 'coord : '
2042       do i=1,inb
2043          !write(2,887) 'bx:',i,bx(i)
2044          !write(2,887) 'by:',i,by(i)
2045          !write(2,887) 'bz:',i,bz(i)
2046       enddo
2047       !write(2,*) 'After sheetforce1'
2048       do i=1,inb
2049          do k=1,12
2050             !write(2,888) 'shefx :',i,k,shefx(i,k)
2051             !write(2,888) 'shefy :',i,k,shefy(i,k)
2052             !write(2,888) 'shefz :',i,k,shefz(i,k)
2053          enddo
2054       enddo
2055
2056       call sheetforce5
2057
2058       !write(2,*) 'After sheetforce5'
2059       do i=1,inb
2060          do k=1,12
2061             !write(2,888) 'shefx :',i,k,shefx(i,k)
2062             !write(2,888) 'shefy :',i,k,shefy(i,k)
2063             !write(2,888) 'shefz :',i,k,shefz(i,k)
2064          enddo
2065       enddo
2066
2067       call sheetforce6
2068
2069       !write(2,*) 'After sheetforce6'
2070       do i=1,inb
2071          do k=1,12
2072             !write(2,888) 'shefx :',i,k,shefx(i,k)
2073             !write(2,888) 'shefy :',i,k,shefy(i,k)
2074             !write(2,888) 'shefz :',i,k,shefz(i,k)
2075          enddo
2076       enddo
2077
2078       call sheetforce11
2079
2080       !write(2,*) 'After sheetforce11'
2081       do i=1,inb
2082          do k=1,12
2083             !write(2,888) 'shefx :',i,k,shefx(i,k)
2084             !write(2,888) 'shefy :',i,k,shefy(i,k)
2085             !write(2,888) 'shefz :',i,k,shefz(i,k)
2086          enddo
2087       enddo
2088
2089       call sheetforce12
2090
2091       !write(2,*) 'After sheetforce12'
2092       do i=1,inb
2093          do k=1,12
2094             !write(2,888) 'shefx :',i,k,shefx(i,k)
2095             !write(2,888) 'shefy :',i,k,shefy(i,k)
2096             !write(2,888) 'shefz :',i,k,shefz(i,k)
2097          enddo
2098       enddo
2099
2100       do i=1,inb
2101          do k=1,12
2102             shetfx(i)=shetfx(i)+shefx(i,k)
2103             shetfy(i)=shetfy(i)+shefy(i,k)
2104             shetfz(i)=shetfz(i)+shefz(i,k)
2105          enddo
2106       enddo
2107       !write(2,*) 'Beta Finished'
2108       do i=1,inb
2109          !write(2,889) 'shetfx : ',i,shetfx(i)
2110          !write(2,889) 'shetfy : ',i,shetfy(i)
2111          !write(2,889) 'shetfz : ',i,shetfz(i)
2112       enddo
2113
2114       return
2115       end
2116 C     end sheetforce
2117 c-------------------------------------------------------------------------------
2118       subroutine sheetene(nca,wshet)
2119       implicit none
2120       integer maxca
2121       parameter(maxca=800)
2122       real*8 dfa_cutoff,dfa_cutoff_delta
2123       parameter(dfa_cutoff=15.5d0)
2124       parameter(dfa_cutoff_delta=0.5d0)
2125 cc******************************************************************************
2126
2127 c      real*8 dfaexp(15001)
2128       real*8 dtmp1, dtmp2, dtmp3
2129
2130       real*8 vbet(maxca,maxca)
2131       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2132       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2133       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2134       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2135       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2136       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2137       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2138       real*8 cph(maxca),cth(maxca)
2139       real*8 rx(maxca,maxca)
2140       real*8 ry(maxca,maxca),rz(maxca,maxca)
2141       real*8 bx(maxca),by(maxca),bz(maxca)
2142       real*8 dis(maxca,maxca)
2143       real*8 ulcos(maxca)
2144 cc**********************************************************************
2145       real*8 astx(maxca),asty(maxca),astz(maxca)
2146       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2147       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2148       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2149       real*8 sth(maxca)
2150       real*8 wshet(maxca,maxca)
2151       real*8 dp45, dm45, w_beta
2152       real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2153       integer nca
2154       integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2155       real*8 uum, uup
2156       real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2157
2158       real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2159       common /shetf/  shetfx,shetfy,shetfz
2160
2161       common /sheca/    bx,by,bz
2162       common /phys1/    inb,nmax,iselect
2163       common /kyori2/   dis
2164       common /difvec/   rx,ry,rz
2165       common /coscos/   cph,cth
2166       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2167      &     c00,s00,ulnex,dnex
2168       common /sheconst/ dp45,dm45,w_beta
2169       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2170       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2171       common /shee/    vbeta,vbet,vbetp,vbetm
2172       common /ulang/    ulcos
2173 cc**********************************************************************
2174       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2175      &     astmmz,astm3x,astm3y,astm3z
2176       common /sinsin/   sth
2177
2178       real*8 r_pair_mat(maxca,maxca)
2179       real*8 e_gcont,fprim_gcont,de_gcont
2180 ci      integer istrand(maxca,maxca)
2181 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2182 ci      common  /shetest/ istrand,istrand_p,istrand_m
2183       common /beta_p/ r_pair_mat
2184 C-------------------------------------------------------------------------------
2185       r_pair_mat = 0.0d0
2186       do i=1,inb
2187          do j=1,inb
2188             r_pair_mat(i,j)=wshet(i,j)
2189 c            write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2190          enddo
2191       enddo
2192 c      stop
2193 c
2194       vbeta=0.0D0
2195       vbetp=0.0D0
2196       vbetm=0.0D0
2197
2198       do i=1,inb-7
2199          do j=i+4,inb-3
2200
2201             if (dis(i,j).lt.dfa_cutoff) then
2202             call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2203      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
2204
2205
2206             ip=i+1
2207             ipp=i+2
2208             jp=j+1
2209             jpp=j+2
2210 cc**********************************************************************
2211             y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2212      &           +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2213             y1=-0.5d0*y1/dca
2214             y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2215      &           +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2216             y2=-0.5d0*y2/dnex
2217
2218 cdebug            y2=0
2219
2220             y=y1+y2
2221
2222 ci           if(y.ge.-4) then
2223 ci              istrand(i,j)=1
2224 ci           else
2225 ci              istrand(i,j)=0
2226 ci           endif
2227
2228 ci           if(istrand(i,j).eq.1) then
2229
2230             yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2231             yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2232
2233
2234             pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2235      $           +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2236             pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2237      $           +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2238             pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2239      $           +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2240             pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2241      $           +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2242
2243            yshe1=pin1(i,j)**2+pin2(i,j)**2
2244            yshe1=-0.5d0*yshe1/dshe
2245            yshe2=pin3(i,j)**2+pin4(i,j)**2
2246            yshe2=-0.5d0*yshe2/dshe
2247
2248 ci              if((yshe1+yshe2).ge.-4) then
2249 ci                 istrand_p(i,j)=1
2250 ci              else
2251 ci                 istrand_p(i,j)=0
2252 ci              endif
2253
2254
2255 C            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2256 C            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2257 C            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2258 C            write(*,*) 'dis(i,j):',i,j,dis(i,j)
2259 C            write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2260 C            write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2261 C            write(*,*) 'pin1:',pin1(i,j)
2262 C            write(*,*) 'pin2:',pin2(i,j)
2263 C            write(*,*) 'pin3:',pin3(i,j)
2264 C            write(*,*) 'pin4:',pin4(i,j)
2265
2266 C            write(*,*) 'y:',y
2267 C            write(*,*) 'yy1:',yy1
2268 C            write(*,*) 'yy2:',yy2
2269 C            write(*,*) 'yshe1:',yshe1
2270 C            write(*,*) 'yshe2:',yshe2
2271 c
2272
2273 ci           if (istrand_p(i,j).eq.1) then
2274
2275 cd           yy1=0
2276 cd           yy2=0
2277 cd           yshe1=0
2278 cd           yshe2=0
2279            dtmp1 = y+yy1+yshe1
2280            dtmp2 = y+yy2+yshe2
2281            dtmp3 = y+yy1+yy2+yshe1+yshe2
2282
2283 C            write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2284 C            write(*,*)'2', y,yy1,yy2
2285 C            write(*,*)'3', yshe1,yshe2
2286
2287 cc           if (dtmp3.le.-35.0d0) then
2288 c              vbetap(i,j)=-dp45*exp(dtmp3)
2289 cc              vbetap(i,j)=0.0d0
2290 cc           else
2291 c              vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2292               vbetap(i,j)=-dp45*exp(dtmp3)
2293 cc           end if
2294
2295 cc           if (dtmp1.le.-35.0d0) then
2296 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2297 cc              vbetap1(i,j)=0.0d0
2298 cc           else
2299 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2300 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2301                vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2302 cc           end if
2303
2304 cc           if (dtmp2.le.-35.0d0) then
2305 C              vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2306 cc              vbetap2(i,j)=0.0d0
2307 cc           else
2308 c              vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2309 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2310               vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2311 cc           end if
2312
2313 c           vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2314 c           vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2315 c           vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2316
2317 !           write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2318 !           write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2319
2320 ci           elseif (istrand_p(i,j).eq.0)then
2321 ci            vbetap(i,j)=0
2322 ci            vbetap1(i,j)=0
2323 ci            vbetap2(i,j)=0
2324 ci           endif
2325
2326            yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2327            yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2328
2329            pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2330      $          +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2331            pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2332      $          +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2333            pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2334      $          +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2335            pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2336      $          +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2337
2338            yshe1=pina1(i,j)**2+pina2(i,j)**2
2339            yshe1=-0.5d0*yshe1/dshe
2340            yshe2=pina3(i,j)**2+pina4(i,j)**2
2341            yshe2=-0.5d0*yshe2/dshe
2342
2343 ci              if((yshe1+yshe2).ge.-4) then
2344 ci                 istrand_m(i,j)=1
2345 ci              else
2346 ci                 istrand_m(i,j)=0
2347 ci              endif
2348
2349
2350 C            write(*,*) 'pina1:',pina1(i,j)
2351 C            write(*,*) 'pina2:',pina2(i,j)
2352 C            write(*,*) 'pina3:',pina3(i,j)
2353 C            write(*,*) 'pina4:',pina4(i,j)
2354 C            write(*,*) 'yshe1:',yshe1
2355 C            write(*,*) 'yshe2:',yshe2
2356 C            write(*,*) 'dshe:',dshe
2357
2358 ci           if (istrand_m(i,j).eq.1) then
2359
2360 cd           yy1=0
2361 cd           yy2=0
2362 cd           yshe1=0
2363 cd           yshe2=0
2364
2365            dtmp3=y+yy1+yy2+yshe1+yshe2
2366            dtmp1=y+yy1+yshe1
2367            dtmp2=y+yy2+yshe2
2368
2369 cc           if(dtmp3 .le. -35.0d0) then
2370 c              vbetam(i,j)=-dm45*exp(dtmp3)
2371 cc              vbetam(i,j)=0.0d0
2372 cc           else
2373 c              vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2374               vbetam(i,j)=-dm45*exp(dtmp3)
2375 cc           end if
2376
2377 cc           if(dtmp1 .le. -35.0d0) then
2378 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2379 cc               vbetam1(i,j)=0.0d0
2380 cc           else
2381 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2382 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2383                vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2384 cc           end if
2385
2386 cc           if(dtmp2.le.-35.0d0) then
2387 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2388 cc              vbetam2(i,j)=0.0d0
2389 cc           else
2390 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2391 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2392               vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2393 cc           end if
2394
2395 ci           elseif (istrand_m(i,j).eq.0)then
2396 ci            vbetam(i,j)=0
2397 ci            vbetam1(i,j)=0
2398 ci            vbetam2(i,j)=0
2399 ci           endif
2400
2401
2402 c           vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2403 c           vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2404 c           vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2405
2406 !           write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2407 !           write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2408
2409            uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2410            uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2411
2412 c           write(*,*) 'uup,uum:', uup, uum
2413
2414 c           uup=vbetap1(i,j)+vbetap2(i,j)
2415 c           uum=vbetam1(i,j)+vbetam2(i,j)
2416
2417            vbet(i,j)=uup+uum
2418            vbetp=vbetp+uup
2419            vbetm=vbetm+uum
2420            vbeta=vbeta+vbet(i,j)*e_gcont
2421
2422
2423            if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
2424 c gradient correction from gcont
2425              de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
2426              shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
2427              shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
2428              shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
2429
2430              shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
2431              shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
2432              shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
2433
2434 c energy correction from gcont
2435              vbet(i,j)=vbet(i,j)*e_gcont
2436              vbetap(i,j)=vbetap(i,j)*e_gcont
2437              vbetap1(i,j)=vbetap1(i,j)*e_gcont
2438              vbetap2(i,j)=vbetap2(i,j)*e_gcont
2439              vbetam(i,j)=vbetam(i,j)*e_gcont
2440              vbetam1(i,j)=vbetam1(i,j)*e_gcont
2441              vbetam2(i,j)=vbetam2(i,j)*e_gcont
2442            endif
2443
2444
2445 ci         elseif(istrand(i,j).eq.0)then
2446 ci           vbet(i,j)=0
2447 ci         endif
2448
2449 c           write(*,*) 'uup,uum:',uup,uum
2450 c           write(*,*) 'vbetap(i,j):',vbetap(i,j)
2451 c           write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2452 c           write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2453 c           write(*,*) 'vbetam(i,j):',vbetam(i,j)
2454 c           write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2455 c           write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2456 c           write(*,*) 'uup:',uup
2457 c           write(*,*) 'uum:',uum
2458 c           write(*,*) 'vbetp:',vbetp
2459 c           write(*,*) 'vbetm:',vbetm
2460 c           write(*,*) 'vbet(i,j):',vbet(i,j)
2461 c           stop
2462
2463             else
2464               vbetap(i,j)=0
2465               vbetap1(i,j)=0
2466               vbetap2(i,j)=0
2467               vbetam(i,j)=0
2468               vbetam1(i,j)=0
2469               vbetam2(i,j)=0
2470               vbet(i,j)=0
2471             endif
2472         enddo
2473       enddo
2474
2475 !      do i=1,inb-7
2476 !         do j=i+4,inb-3
2477 !            write(*,*) 'I,J:', i,j
2478 !            write(*,*) 'vbetap(i,j):',vbetap(i,j)
2479 !            write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2480 !            write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2481 !            write(*,*) 'vbetam(i,j):',vbetam(i,j)
2482 !            write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2483 !            write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2484 !            write(*,*) 'vbet(i,j):',vbet(i,j)
2485 !         enddo
2486 !      enddo
2487
2488       return
2489       end
2490 c-------------------------------------------------------------------------------
2491       subroutine sheetforce1
2492       implicit none
2493       integer maxca
2494       parameter(maxca=800)
2495       real*8 dfa_cutoff,dfa_cutoff_delta
2496       parameter(dfa_cutoff=15.5d0)
2497       parameter(dfa_cutoff_delta=0.5d0)
2498 cc**********************************************************************
2499       real*8 vbet(maxca,maxca)
2500       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2501       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2502       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2503       real*8 cph(maxca),cth(maxca)
2504       real*8 rx(maxca,maxca)
2505       real*8 ry(maxca,maxca),rz(maxca,maxca)
2506       real*8 bx(maxca),by(maxca),bz(maxca)
2507       real*8 dis(maxca,maxca)
2508       real*8 shefx(maxca,12)
2509       real*8 shefy(maxca,12),shefz(maxca,12)
2510       real*8 atx(maxca),aty(maxca),atz(maxca)
2511       real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2512       real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2513       real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2514       real*8 apx(maxca),apy(maxca),apz(maxca)
2515       real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2516       real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2517       real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2518       real*8 ulcos(maxca)
2519       real*8 astx(maxca),asty(maxca),astz(maxca)
2520       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2521       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2522       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2523       real*8 sth(maxca)
2524       real*8 w_beta,dp45, dm45
2525       real*8 vbeta, vbetp, vbetm
2526       real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2527      $     c00,s00,ulnex,dnex
2528       integer inb,nmax,iselect
2529
2530       common /phys1/     inb,nmax,iselect
2531       common /kyori2/    dis
2532       common /difvec/   rx,ry,rz
2533       common /coscos/   cph,cth
2534       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2535      $     c00,s00,ulnex,dnex
2536       common /sheconst/ dp45,dm45,w_beta
2537       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2538       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2539      $     atmmz,atm3x,atm3y,atm3z
2540       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2541      $     apmmz,apm3x,apm3y,apm3z
2542       common /shef/   shefx,shefy,shefz
2543       common /shee/   vbeta,vbet,vbetp,vbetm
2544       common /ulang/    ulcos
2545 c     c**********************************************************************
2546       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2547      $     astmmz,astm3x,astm3y,astm3z
2548       common /sinsin/   sth
2549 C--------------------------------------------------------------------------------
2550 c     local variables
2551       integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2552       real*8  c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2553       real*8  c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2554       real*8  c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2555       real*8  dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2556       real*8 e_gcont,fprim_gcont
2557 C--------------------------------------------------------------------------------
2558       do i=4,inb-4
2559          im3=i-3
2560          imm=i-2
2561          im=i-1
2562          c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2563          v1=0.0D0
2564          do j=i+1,inb-3
2565             v1=v1+vbet(im3,j)
2566          enddo
2567          cc1=(ulcos(imm)-ulnex)/dnex
2568          dmm=cc1/(dis(imm,im)*dis(im,i))
2569          dmm__=cc1*ulcos(imm)/dis(im,i)**2
2570          fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2571          fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2572          fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2573 cd         fx=0
2574 cd         fy=0
2575 cd         fz=0
2576          fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2577          fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2578          fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2579          shefx(i,1)=fx*v1
2580          shefy(i,1)=fy*v1
2581          shefz(i,1)=fz*v1
2582       enddo
2583
2584       do i=3,inb-5
2585          imm=i-2
2586          im=i-1
2587          ip=i+1
2588          c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2589          v2=0.0D0
2590          do j=i+2,inb-3
2591             v2=v2+vbet(imm,j)
2592          enddo
2593          cc1=(ulcos(imm)-ulnex)/dnex
2594          cc2=(ulcos(im)-ulnex)/dnex
2595          dmm1=cc1/(dis(imm,im)*dis(im,i))
2596          dmm2=cc2/(dis(im,i)*dis(i,ip))
2597          dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2598          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2599          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2600 cc**********************************************************************
2601          fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2602      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2603          fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2604      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2605          fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2606      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2607 cd         fx=0
2608 cd         fy=0
2609 cd         fz=0
2610          fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2611          fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2612          fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2613          shefx(i,2)=fx*v2
2614          shefy(i,2)=fy*v2
2615          shefz(i,2)=fz*v2
2616       enddo
2617       do i=2,inb-6
2618          im=i-1
2619          ip=i+1
2620          ipp=i+2
2621          c3=(cth(im)*c00+sth(im)*s00-1)/dca
2622          v3=0.0D0
2623          do j=i+3,inb-3
2624             v3=v3+vbet(im,j)
2625          enddo
2626          cc2=(ulcos(im)-ulnex)/dnex
2627          cc3=(ulcos(i)-ulnex)/dnex
2628          dmm2=cc2/(dis(im,i)*dis(i,ip))
2629          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2630          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2631          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2632          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2633          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2634      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2635          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2636      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2637          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2638      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2639 cd         fx=0
2640 cd         fy=0
2641 cd         fz=0
2642          fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2643          fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2644          fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2645          shefx(i,3)=fx*v3
2646          shefy(i,3)=fy*v3
2647          shefz(i,3)=fz*v3
2648       enddo
2649       do i=1,inb-7
2650          ip=i+1
2651          ipp=i+2
2652          c4=(cth(i)*c00+sth(i)*s00-1)/dca
2653          v4=0.0D0
2654          do j=i+4,inb-3
2655             v4=v4+vbet(i,j)
2656          enddo
2657          cc3=(ulcos(i)-ulnex)/dnex
2658          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2659          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2660          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2661          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2662          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2663 cd         fx=0
2664 cd         fy=0
2665 cd         fz=0
2666          fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2667          fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2668          fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2669          shefx(i,4)=fx*v4
2670          shefy(i,4)=fy*v4
2671          shefz(i,4)=fz*v4
2672       enddo
2673       do j=8,inb
2674          jm3=j-3
2675          jmm=j-2
2676          jm=j-1
2677          c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2678          v7=0.0D0
2679          do i=1,j-7
2680             v7=v7+vbet(i,jm3)
2681          enddo
2682          cc7=(ulcos(jmm)-ulnex)/dnex
2683          dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2684          dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2685          fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2686          fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2687          fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2688 cd         fx=0
2689 cd         fy=0
2690 cd         fz=0
2691          fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2692          fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2693          fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2694          shefx(j,7)=fx*v7
2695          shefy(j,7)=fy*v7
2696          shefz(j,7)=fz*v7
2697       enddo
2698       do j=7,inb-1
2699          jm=j-1
2700          jmm=j-2
2701          jp=j+1
2702          c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2703          v8=0.0D0
2704          do i=1,j-6
2705             v8=v8+vbet(i,jmm)
2706          enddo
2707          cc7=(ulcos(jmm)-ulnex)/dnex
2708          cc8=(ulcos(jm)-ulnex)/dnex
2709          dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2710          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2711          dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2712          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2713          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2714          fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2715      $        -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2716          fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2717      $        -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2718          fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2719      $        -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2720 cd         fx=0
2721 cd         fy=0
2722 cd         fz=0
2723          fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2724          fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2725          fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2726          shefx(j,8)=fx*v8
2727          shefy(j,8)=fy*v8
2728          shefz(j,8)=fz*v8
2729       enddo
2730
2731       do j=6,inb-2
2732          jm=j-1
2733          jp=j+1
2734          jpp=j+2
2735          c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2736          v9=0.0D0
2737          do i=1,j-5
2738             v9=v9+vbet(i,jm)
2739          enddo
2740          cc8=(ulcos(jm)-ulnex)/dnex
2741          cc9=(ulcos(j)-ulnex)/dnex
2742          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2743          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2744          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2745          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2746          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2747          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2748      $        -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2749          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2750      $        -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2751          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2752      $        -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2753 cd         fx=0
2754 cd         fy=0
2755 cd         fz=0
2756          fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2757          fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2758          fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2759          shefx(j,9)=fx*v9
2760          shefy(j,9)=fy*v9
2761          shefz(j,9)=fz*v9
2762       enddo
2763
2764       do j=5,inb-3
2765          jp=j+1
2766          jpp=j+2
2767          c10=(cth(j)*c00+sth(j)*s00-1)/dca
2768          v10=0.0D0
2769          do i=1,j-4
2770             v10=v10+vbet(i,j)
2771          enddo
2772          cc9=(ulcos(j)-ulnex)/dnex
2773          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2774          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2775          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2776          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2777          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2778 cd         fx=0
2779 cd         fy=0
2780 cd         fz=0
2781          fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2782          fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2783          fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2784          shefx(j,10)=fx*v10
2785          shefy(j,10)=fy*v10
2786          shefz(j,10)=fz*v10
2787       enddo
2788
2789       return
2790       end
2791 c----------------------------------------------------------------------------
2792       subroutine sheetforce5
2793       implicit none
2794       integer maxca
2795       parameter(maxca=800)
2796       real*8 dfa_cutoff,dfa_cutoff_delta
2797       parameter(dfa_cutoff=15.5d0)
2798       parameter(dfa_cutoff_delta=0.5d0)
2799 cc**********************************************************************
2800       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2801       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2802       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2803       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2804       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2805       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2806       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2807       real*8 rx(maxca,maxca)
2808       real*8 ry(maxca,maxca),rz(maxca,maxca)
2809       real*8 bx(maxca),by(maxca),bz(maxca)
2810       real*8 dis(maxca,maxca)
2811       real*8 shefx(maxca,12),shefy(maxca,12)
2812       real*8 shefz(maxca,12)
2813       real*8 dp45,dm45,w_beta
2814       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2815      $     c00,s00,ulnex,dnex
2816       integer    inb,nmax,iselect
2817 cc**********************************************************************
2818       common /phys1/     inb,nmax,iselect
2819       common /kyori2/    dis
2820       common /difvec/   rx,ry,rz
2821       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2822      $     c00,s00,ulnex,dnex
2823       common /sheconst/ dp45,dm45,w_beta
2824       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2825       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2826       common /shef/   shefx,shefy,shefz
2827 ci      integer istrand(maxca,maxca)
2828 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2829 ci      common  /shetest/ istrand,istrand_p,istrand_m
2830 c********************************************************************************
2831 c     local variables
2832       integer i,imm,im,jp,jpp,j
2833       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2834       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2835       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2836       real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2837       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2838       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2839       real*8 e_gcont,fprim_gcont
2840 c********************************************************************************
2841       do i=3,inb-5
2842          imm=i-2
2843          im=i-1
2844          do j=i+2,inb-3
2845
2846             if (dis(imm,j).lt.dfa_cutoff) then
2847             call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2848      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
2849
2850             jp=j+1
2851             jpp=j+2
2852
2853 ci            if(istrand(imm,j).eq.1
2854 ci     &   .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2855
2856
2857             yy1=-(dis(i,jpp)-ulhb)/dlhb
2858             y1x=rx(jpp,i)/dis(i,jpp)
2859             y1y=ry(jpp,i)/dis(i,jpp)
2860             y1z=rz(jpp,i)/dis(i,jpp)
2861             y11x=yy1*y1x
2862             y11y=yy1*y1y
2863             y11z=yy1*y1z
2864
2865             yy33=1.0D0/(dis(im,jp)*dis(im,i))
2866             yyy3=pin1(imm,j)/(dis(im,i)**2)
2867             yy3=-pin1(imm,j)/dshe
2868             y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2869             y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2870             y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2871
2872             yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2873             yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2874             yyy4b=pin3(imm,j)/(dis(im,i)**2)
2875             yy4=-pin3(imm,j)/dshe
2876             y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2877      $           -yyy4b*rx(im,i))*yy4
2878             y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2879      $           -yyy4b*ry(im,i))*yy4
2880             y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2881      $           -yyy4b*rz(im,i))*yy4
2882
2883
2884             yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2885             yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2886             yy5=-pin4(imm,j)/dshe
2887             y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2888             y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2889             y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2890
2891             sx=y11x+y3x+y4x+y5x
2892             sy=y11y+y3y+y4y+y5y
2893             sz=y11z+y3z+y4z+y5z
2894
2895             sx1=y3x
2896             sy1=y3y
2897             sz1=y3z
2898             sx2=y11x+y4x+y5x
2899             sy2=y11y+y4y+y5y
2900             sz2=y11z+y4z+y5z
2901
2902             shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2903      $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2904             shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2905      $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2906             shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2907      $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2908
2909 !            shefx(i,5)=shefx(i,5)
2910 !     $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2911 !            shefy(i,5)=shefy(i,5)
2912 !     $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2913 !            shefz(i,5)=shefz(i,5)
2914 !     $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2915
2916             yy6=-(dis(i,jp)-uldhb)/dldhb
2917             y6x=rx(jp,i)/dis(i,jp)
2918             y6y=ry(jp,i)/dis(i,jp)
2919             y6z=rz(jp,i)/dis(i,jp)
2920             y66x=yy6*y6x
2921             y66y=yy6*y6y
2922             y66z=yy6*y6z
2923
2924             yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2925             yyy8=pina1(imm,j)/(dis(im,i)**2)
2926             yy8=-pina1(imm,j)/dshe
2927             y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2928             y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2929             y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2930
2931             yy99=1.0D0/(dis(jp,i)*dis(im,i))
2932             yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2933             yyy9b=pina3(imm,j)/(dis(im,i)**2)
2934             yy9=-pina3(imm,j)/dshe
2935             y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2936      $           -yyy9b*rx(im,i))*yy9
2937             y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2938      $           -yyy9b*ry(im,i))*yy9
2939             y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2940      $           -yyy9b*rz(im,i))*yy9
2941
2942             yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2943             yyy10=pina4(imm,j)/(dis(jp,i)**2)
2944             yy10=-pina4(imm,j)/dshe
2945             y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2946             y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2947             y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2948
2949             sx=y66x+y8x+y9x+y10x
2950             sy=y66y+y8y+y9y+y10y
2951             sz=y66z+y8z+y9z+y10z
2952
2953             sx1=y8x
2954             sy1=y8y
2955             sz1=y8z
2956             sx2=y66x+y9x+y10x
2957             sy2=y66y+y9y+y10y
2958             sz2=y66z+y9z+y10z
2959
2960             shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
2961      $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2962            shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
2963      $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2964             shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
2965      $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2966
2967 !            shefx(i,5)=shefx(i,5)
2968 !     $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2969 !            shefy(i,5)=shefy(i,5)
2970 !     $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2971 !            shefz(i,5)=shefz(i,5)
2972 !     $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2973
2974              endif
2975 ci          endif
2976
2977          enddo
2978       enddo
2979
2980       return
2981       end
2982 c--------------------------------------------------------------------------c
2983       subroutine sheetforce6
2984       implicit none
2985       integer maxca
2986       parameter(maxca=800)
2987       real*8 dfa_cutoff,dfa_cutoff_delta
2988       parameter(dfa_cutoff=15.5d0)
2989       parameter(dfa_cutoff_delta=0.5d0)
2990 cc**********************************************************************
2991       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2992       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2993       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2994       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2995       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2996       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2997       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2998       real*8 rx(maxca,maxca)
2999       real*8 ry(maxca,maxca),rz(maxca,maxca)
3000       real*8 bx(maxca),by(maxca),bz(maxca)
3001       real*8 dis(maxca,maxca)
3002       real*8 shefx(maxca,12),shefy(maxca,12)
3003       real*8 shefz(maxca,12)
3004       real*8 dp45,dm45,w_beta
3005       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3006      $     c00,s00,ulnex,dnex
3007       integer    inb,nmax,iselect
3008 cc**********************************************************************
3009       common /phys1/     inb,nmax,iselect
3010       common /kyori2/    dis
3011       common /difvec/   rx,ry,rz
3012       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3013      $     c00,s00,ulnex,dnex
3014       common /sheconst/ dp45,dm45,w_beta
3015       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3016       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3017       common /shef/   shefx,shefy,shefz
3018 ci      integer istrand(maxca,maxca)
3019 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3020 ci      common  /shetest/ istrand,istrand_p,istrand_m
3021 cc**********************************************************************
3022 C     local variables
3023       integer  i,imm,im,jp,jpp,j,ip
3024       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3025       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3026       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3027       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3028       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
3029       real*8  yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
3030       real*8 e_gcont,fprim_gcont
3031 C********************************************************************************
3032       do i=2,inb-6
3033          ip=i+1
3034          im=i-1
3035          do j=i+3,inb-3
3036
3037             if (dis(im,j).lt.dfa_cutoff) then
3038             call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3039      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3040
3041             jp=j+1
3042             jpp=j+2
3043
3044 ci        if(istrand(im,j).eq.1
3045 ci     &    .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
3046
3047
3048             yy1=-(dis(i,jp)-ulhb)/dlhb
3049             y1x=rx(jp,i)/dis(i,jp)
3050             y1y=ry(jp,i)/dis(i,jp)
3051             y1z=rz(jp,i)/dis(i,jp)
3052             y11x=yy1*y1x
3053             y11y=yy1*y1y
3054             y11z=yy1*y1z
3055
3056             yy33=1.0D0/(dis(i,jp)*dis(i,ip))
3057             yyy3a=pin1(im,j)/(dis(i,jp)**2)
3058             yyy3b=pin1(im,j)/(dis(i,ip)**2)
3059             yy3=-pin1(im,j)/dshe
3060             y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
3061      $           +yyy3b*rx(i,ip))*yy3
3062             y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
3063      $           +yyy3b*ry(i,ip))*yy3
3064             y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
3065      $           +yyy3b*rz(i,ip))*yy3
3066
3067             yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3068             yyy4=pin2(im,j)/(dis(i,jp)**2)
3069             yy4=-pin2(im,j)/dshe
3070             y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3071             y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3072             y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3073
3074             yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3075             yyy5=pin3(im,j)/(dis(i,ip)**2)
3076             yy5=-pin3(im,j)/dshe
3077             y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3078             y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3079             y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3080
3081             sx=y11x+y3x+y4x+y5x
3082             sy=y11y+y3y+y4y+y5y
3083             sz=y11z+y3z+y4z+y5z
3084
3085             sx1=y11x+y3x+y4x
3086             sy1=y11y+y3y+y4y
3087             sz1=y11z+y3z+y4z
3088             sx2=y5x
3089             sy2=y5y
3090             sz2=y5z
3091
3092             shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3093      $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3094             shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3095      $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3096             shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3097      $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3098 !            shefx(i,6)=shefx(i,6)
3099 !     $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3100 !            shefy(i,6)=shefy(i,6)
3101 !     $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3102 !            shefz(i,6)=shefz(i,6)
3103 !     $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3104
3105             yy6=-(dis(jpp,i)-uldhb)/dldhb
3106             y6x=rx(jpp,i)/dis(jpp,i)
3107             y6y=ry(jpp,i)/dis(jpp,i)
3108             y6z=rz(jpp,i)/dis(jpp,i)
3109             y66x=yy6*y6x
3110             y66y=yy6*y6y
3111             y66z=yy6*y6z
3112
3113             yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3114             yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3115             yyy8b=pina1(im,j)/(dis(i,ip)**2)
3116             yy8=-pina1(im,j)/dshe
3117             y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3118      $           +yyy8b*rx(i,ip))*yy8
3119             y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3120      $           +yyy8b*ry(i,ip))*yy8
3121             y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3122      $           +yyy8b*rz(i,ip))*yy8
3123
3124             yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3125             yyy9=pina2(im,j)/(dis(i,jpp)**2)
3126             yy9=-pina2(im,j)/dshe
3127             y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3128             y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3129             y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3130
3131             yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3132             yyy10=pina3(im,j)/(dis(i,ip)**2)
3133             yy10=-pina3(im,j)/dshe
3134             y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3135             y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3136             y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3137
3138             sx=y66x+y8x+y9x+y10x
3139             sy=y66y+y8y+y9y+y10y
3140             sz=y66z+y8z+y9z+y10z
3141
3142             sx1=y66x+y8x+y9x
3143             sy1=y66y+y8y+y9y
3144             sz1=y66z+y8z+y9z
3145             sx2=y10x
3146             sy2=y10y
3147             sz2=y10z
3148
3149             shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3150      $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3151            shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3152      $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3153             shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3154      $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3155
3156 !            shefx(i,6)=shefx(i,6)
3157 !     $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3158 !           shefy(i,6)=shefy(i,6)
3159 !     $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3160 !            shefz(i,6)=shefz(i,6)
3161 !     $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3162
3163             endif
3164 ci         endif
3165
3166          enddo
3167       enddo
3168
3169       return
3170       end
3171 c-----------------------------------------------------------------------
3172       subroutine sheetforce11
3173       implicit none
3174       integer maxca
3175       parameter(maxca=800)
3176       real*8 dfa_cutoff,dfa_cutoff_delta
3177       parameter(dfa_cutoff=15.5d0)
3178       parameter(dfa_cutoff_delta=0.5d0)
3179 cc**********************************************************************
3180       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3181       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3182       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3183       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3184       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3185       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3186       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3187       real*8 rx(maxca,maxca)
3188       real*8 ry(maxca,maxca),rz(maxca,maxca)
3189       real*8 bx(maxca),by(maxca),bz(maxca)
3190       real*8 dis(maxca,maxca)
3191       real*8 shefx(maxca,12),shefy(maxca,12)
3192       real*8 shefz(maxca,12)
3193       real*8 dp45,dm45,w_beta
3194       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3195      $     c00,s00,ulnex,dnex
3196       integer    inb,nmax,iselect
3197 cc**********************************************************************
3198       common /phys1/     inb,nmax,iselect
3199       common /kyori2/    dis
3200       common /difvec/   rx,ry,rz
3201       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3202      $     c00,s00,ulnex,dnex
3203       common /sheconst/ dp45,dm45,w_beta
3204       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3205       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3206       common /shef/   shefx,shefy,shefz
3207 ci      integer istrand(maxca,maxca)
3208 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3209 ci      common  /shetest/ istrand,istrand_p,istrand_m
3210 C********************************************************************************
3211 C     local variables
3212       integer  j,jm,jmm,ip,i,ipp
3213       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3214       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3215       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3216       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3217       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3218       real*8  yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3219       real*8 e_gcont,fprim_gcont
3220 C********************************************************************************
3221
3222       do j=7,inb-1
3223          jm=j-1
3224          jmm=j-2
3225          do i=1,j-6
3226
3227             if (dis(i,jmm).lt.dfa_cutoff) then
3228             call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3229      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3230
3231             ip=i+1
3232             ipp=i+2
3233
3234 ci            if(istrand(i,jmm).eq.1
3235 ci     &   .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3236
3237
3238             yy1=-(dis(ipp,j)-ulhb)/dlhb
3239             y1x=rx(ipp,j)/dis(ipp,j)
3240             y1y=ry(ipp,j)/dis(ipp,j)
3241             y1z=rz(ipp,j)/dis(ipp,j)
3242             y11x=yy1*y1x
3243             y11y=yy1*y1y
3244             y11z=yy1*y1z
3245
3246             yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3247             yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3248             yy3=-pin2(i,jmm)/dshe
3249             y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3250             y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3251             y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3252
3253             yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3254             yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3255             yy4=-pin3(i,jmm)/dshe
3256             y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3257             y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3258             y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3259
3260             yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3261             yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3262             yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3263             yy5=-pin4(i,jmm)/dshe
3264             y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3265      $           -yyy5b*rx(jm,j))*yy5
3266             y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3267      $           -yyy5b*ry(jm,j))*yy5
3268             y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3269      $           -yyy5b*rz(jm,j))*yy5
3270
3271             sx=y11x+y3x+y4x+y5x
3272             sy=y11y+y3y+y4y+y5y
3273             sz=y11z+y3z+y4z+y5z
3274
3275             sx1=y3x
3276             sy1=y3y
3277             sz1=y3z
3278             sx2=y11x+y4x+y5x
3279             sy2=y11y+y4y+y5y
3280             sz2=y11z+y4z+y5z
3281
3282             shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3283      $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3284             shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3285      $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3286             shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3287      $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3288
3289 !            shefx(j,11)=shefx(j,11)
3290 !     $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3291 !            shefy(j,11)=shefy(j,11)
3292 !     $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3293 !            shefz(j,11)=shefz(j,11)
3294 !     $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3295
3296             yy6=-(dis(ip,j)-uldhb)/dldhb
3297             y6x=rx(ip,j)/dis(ip,j)
3298             y6y=ry(ip,j)/dis(ip,j)
3299             y6z=rz(ip,j)/dis(ip,j)
3300             y66x=yy6*y6x
3301             y66y=yy6*y6y
3302             y66z=yy6*y6z
3303
3304             yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3305             yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3306             yy8=-pina1(i,jmm)/dshe
3307             y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3308             y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3309             y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3310
3311             yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3312             yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3313             yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3314             yy9=-pina2(i,jmm)/dshe
3315             y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3316      $           -yyy9b*rx(jm,j))*yy9
3317             y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3318      $           -yyy9b*ry(jm,j))*yy9
3319             y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3320      $           -yyy9b*rz(jm,j))*yy9
3321
3322             yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3323             yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3324             yy10=-pina4(i,jmm)/dshe
3325             y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3326             y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3327             y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3328
3329             sx=y66x+y8x+y9x+y10x
3330             sy=y66y+y8y+y9y+y10y
3331             sz=y66z+y8z+y9z+y10z
3332
3333             sx1=y66x+y8x+y9x
3334             sy1=y66y+y8y+y9y
3335             sz1=y66z+y8z+y9z
3336             sx2=y10x
3337             sy2=y10y
3338             sz2=y10z
3339
3340             shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3341      $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3342            shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3343      $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3344             shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3345      $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3346
3347 !            shefx(j,11)=shefx(j,11)
3348 !     $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3349 !            shefy(j,11)=shefy(j,11)
3350 !     $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3351 !            shefz(j,11)=shefz(j,11)
3352 !     $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3353
3354             endif
3355 ci         endif
3356
3357          enddo
3358       enddo
3359
3360       return
3361       end
3362 c-----------------------------------------------------------------------
3363       subroutine sheetforce12
3364       implicit none
3365       integer maxca
3366       parameter(maxca=800)
3367       real*8 dfa_cutoff,dfa_cutoff_delta
3368       parameter(dfa_cutoff=15.5d0)
3369       parameter(dfa_cutoff_delta=0.5d0)
3370 cc**********************************************************************
3371       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3372       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3373       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3374       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3375       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3376       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3377       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3378       real*8 rx(maxca,maxca)
3379       real*8 ry(maxca,maxca),rz(maxca,maxca)
3380       real*8 bx(maxca),by(maxca),bz(maxca)
3381       real*8 dis(maxca,maxca)
3382       real*8 shefx(maxca,12),shefy(maxca,12)
3383       real*8 shefz(maxca,12)
3384       real*8 dp45,dm45,w_beta
3385       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3386      $     c00,s00,ulnex,dnex
3387       integer    inb,nmax,iselect
3388 cc**********************************************************************
3389       common /phys1/     inb,nmax,iselect
3390       common /kyori2/    dis
3391       common /difvec/   rx,ry,rz
3392       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3393      $     c00,s00,ulnex,dnex
3394       common /sheconst/ dp45,dm45,w_beta
3395       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3396       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3397       common /shef/   shefx,shefy,shefz
3398 ci      integer istrand(maxca,maxca)
3399 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3400 ci      common  /shetest/ istrand,istrand_p,istrand_m
3401 cc**********************************************************************
3402 C     local variables
3403       integer j,jm,jmm,ip,i,ipp,jp
3404       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3405       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3406       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3407       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3408       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3409       real*8 e_gcont,fprim_gcont
3410 !c*************************************************************************c
3411       do j=6,inb-2
3412          jp=j+1
3413          jm=j-1
3414          do i=1,j-5
3415
3416             if (dis(i,jm).lt.dfa_cutoff) then
3417             call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3418      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3419
3420             ip=i+1
3421             ipp=i+2
3422
3423 ci            if(istrand(i,jm).eq.1
3424 ci     &   .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3425
3426
3427             yy1=-(dis(ip,j)-ulhb)/dlhb
3428             y1x=rx(ip,j)/dis(ip,j)
3429             y1y=ry(ip,j)/dis(ip,j)
3430             y1z=rz(ip,j)/dis(ip,j)
3431             y11x=y1x*yy1
3432             y11y=y1y*yy1
3433             y11z=y1z*yy1
3434
3435             yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3436             yyy3=pin1(i,jm)/(dis(ip,j)**2)
3437             yy3=-pin1(i,jm)/dshe
3438             y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3439             y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3440             y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3441             yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3442
3443             yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3444             yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3445             yy4=-pin2(i,jm)/dshe
3446             y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3447      $           +yyy4b*rx(j,jp))*yy4
3448             y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3449      $           +yyy4b*ry(j,jp))*yy4
3450             y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3451      $           +yyy4b*rz(j,jp))*yy4
3452
3453             yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3454             yyy5=pin4(i,jm)/(dis(j,jp)**2)
3455             yy5=-pin4(i,jm)/dshe
3456             y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3457             y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3458             y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3459
3460             sx=y11x+y3x+y4x+y5x
3461             sy=y11y+y3y+y4y+y5y
3462             sz=y11z+y3z+y4z+y5z
3463
3464             sx1=y11x+y3x+y4x
3465             sy1=y11y+y3y+y4y
3466             sz1=y11z+y3z+y4z
3467             sx2=y5x
3468             sy2=y5y
3469             sz2=y5z
3470
3471             shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3472      $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3473             shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3474      $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3475             shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3476      $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3477
3478 !            shefx(j,12)=shefx(j,12)
3479 !     $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3480 !            shefy(j,12)=shefy(j,12)
3481 !     $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3482 !            shefz(j,12)=shefz(j,12)
3483 !     $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3484
3485             yy6=-(dis(ipp,j)-uldhb)/dldhb
3486             y6x=rx(ipp,j)/dis(ipp,j)
3487             y6y=ry(ipp,j)/dis(ipp,j)
3488             y6z=rz(ipp,j)/dis(ipp,j)
3489             y66x=yy6*y6x
3490             y66y=yy6*y6y
3491             y66z=yy6*y6z
3492
3493             yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3494             yyy8=pina2(i,jm)/(dis(j,jp)**2)
3495             yy8=-pina2(i,jm)/dshe
3496             y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3497             y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3498             y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3499
3500             yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3501             yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3502             yy9=-pina3(i,jm)/dshe
3503             y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3504             y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3505             y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3506
3507             yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3508             yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3509             yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3510             yy10=-pina4(i,jm)/dshe
3511             y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3512      $           +yyy10b*rx(j,jp))*yy10
3513             y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3514      $           +yyy10b*ry(j,jp))*yy10
3515             y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3516      $           +yyy10b*rz(j,jp))*yy10
3517
3518             sx=y66x+y8x+y9x+y10x
3519             sy=y66y+y8y+y9y+y10y
3520             sz=y66z+y8z+y9z+y10z
3521
3522             sx1=y8x
3523             sy1=y8y
3524             sz1=y8z
3525             sx2=y66x+y9x+y10x
3526             sy2=y66y+y9y+y10y
3527             sz2=y66z+y9z+y10z
3528
3529             shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3530      $           -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3531            shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3532      $           -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3533             shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3534      $           -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3535
3536             endif
3537
3538 ci         endif
3539
3540          ENDDO
3541       ENDDO
3542
3543       RETURN
3544       END
3545 C===============================================================================