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