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