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