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