1 subroutine init_dfa_vars
4 include 'COMMON.INTERACT'
28 CC WEIGHTS for each min
45 C weights for each E term
46 C these should be identical with
53 C precalculate exp table!
56 c dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
62 print *,'ishiftca=',ishiftca,'ilastca=',ilastca
68 subroutine read_dfa_info
70 C read fragment informations
72 implicit real*8 (a-h,o-z)
76 include 'COMMON.SETUP'
79 include 'COMMON.IOUNITS'
80 include 'COMMON.CHAIN'
82 include 'COMMON.FFIELD'
83 include 'COMMON.CONTROL'
86 C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
87 C THIS SHOULD BE MODIFIED!!
94 integer ica1, ica2,ica3,ica4,ica5
95 integer ishell, inca, itmp,iitmp
97 logical lprn /.false./
101 open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
103 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
104 print *,'Error opening dist_dfa.dat file'
107 write(iout,'(a)') 'dist_dfa.dat is opened!'
109 read(iodfa, '(a)') buffer
110 C read number of restraints
111 read(iodfa, *) IDFADIS
112 read(iodfa, *) dis_inc
114 read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
133 C READ ANGLE RESTRAINTS
135 open(iodfa, file='phi_dfa.dat',status='old',err=35)
137 35 write(iout,'(a)') 'Error opening phi_dfa.dat file'
138 print *, 'Error opening phi_dfa.dat file'
142 write(iout,'(a)') 'phi_dfa.dat is opened!'
145 read(iodfa, '(a)') buffer
146 C READ NUMBER OF RESTRAINTS
147 READ(iodfa, *) IDFAPHI
148 read(iodfa,*) phi_inc
150 read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
161 read(iodfa,*) tmp1,tmp2
175 open(iodfa, file='theta_dfa.dat',status='old',err=41)
177 41 write(iout,'(a)') 'Error opening theta_dfa.dat file'
178 print *,'Error opening theta_dfa.dat file'
181 write(iout,'(a)') 'theta_dfa.dat is opened!'
183 read(iodfa, '(a)') buffer
184 C READ NUMBER OF RESTRAINTS
185 READ(iodfa, *) IDFATHE
186 read(iodfa,*) the_inc
189 read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
200 read(iodfa,*) tmp1,tmp2
212 C END of READING ANGLE RESTRAINT!
214 C NUMBER OF NEIGHBOR CAs
215 open(iodfa,file='nei_dfa.dat',status='old',err=37)
217 37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
218 print *,'Error opening nei_dfa.dat file'
221 write(iout,'(a)') 'nei_dfa.dat is opened!'
223 read(iodfa, '(a)') buffer
224 C READ NUMBER OF RESTRAINTS
225 READ(iodfa, *) idfanei
226 read(iodfa,*) nei_inc
229 read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
238 C write(*,*) 'READ NEI:',i,j,fnei(i,j)
248 C END OF NEIGHBORING CA
251 C BETA is not parallel !
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"
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
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 if (nfgtasks.gt.1) then
274 do i=0,max_fg_procs-1
275 idfadis_start_all(j)=0
277 idfaphi_start_all(j)=0
279 idfathe_start_all(j)=0
281 idfanei_start_all(j)=0
284 call MPI_Allgather(idfadis_start,1,MPI_INTEGER,
285 & idfadis_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
286 call MPI_Allgather(idfadis_end,1,MPI_INTEGER,
287 & idfadis_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
288 call MPI_Allgather(idfaphi_start,1,MPI_INTEGER,
289 & idfaphi_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
290 call MPI_Allgather(idfaphi_end,1,MPI_INTEGER,
291 & idfaphi_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
292 call MPI_Allgather(idfathe_start,1,MPI_INTEGER,
293 & idfathe_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
294 call MPI_Allgather(idfathe_end,1,MPI_INTEGER,
295 & idfathe_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
296 call MPI_Allgather(idfanei_start,1,MPI_INTEGER,
297 & idfanei_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
298 call MPI_Allgather(idfanei_end,1,MPI_INTEGER,
299 & idfanei_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
300 if (me.eq.0 .or. out1file) then
301 write (iout,*) "Partitioning of DFA work"
302 write (iout,'(5a10)') 'Rank','DFA_dis','DFA_phi','DFA_the',
305 write (iout,'(i10,8i5)') i,idfadis_start_all(i),
306 & idfadis_end_all(i),idfaphi_start_all(i),
307 & idfaphi_end_all(i),idfathe_start_all(i),
308 & idfathe_end_all(i),idfanei_start_all(i),
325 C READ BETA RESTRAINT
326 if (wdfa_beta.eq.0.0) return
327 open(iodfa, file='beta_dfa.dat',status='old',err=39)
329 39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
330 print *,'Error opening beta_dfa.dat file'
333 write(iout,'(a)') 'beta_dfa.dat is opened!'
335 read(iodfa,'(a)') buffer
337 read(iodfa,*) beta_inc
340 read(iodfa,*) ica1, iitmp
344 c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
349 C END OF BETA RESTRAINT
355 subroutine edfad(edfadis)
357 implicit real*8 (a-h,o-z)
359 include 'COMMON.CHAIN'
360 include 'COMMON.DERIV'
363 double precision edfadis
364 integer i, iatm1, iatm2,idiff
365 double precision ckk, sckk,dist,texp
366 double precision jix,jiy,jiz,ep,fp,scc
371 c write (2,*) "edfad",idfadis_start,idfadis_end
372 do i=idfadis_start,idfadis_end
374 iatm1=idislis(1,i)+ishiftca
375 iatm2=idislis(2,i)+ishiftca
376 idiff = abs(iatm1-iatm2)
378 JIX=c(1,iatm2)-c(1,iatm1)
379 JIY=c(2,iatm2)-c(2,iatm1)
380 JIZ=c(3,iatm2)-c(3,iatm1)
381 DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
394 if (dtmp.ge.15.0d0) then
397 c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
398 texp = exp(-dtmp)/sckk
401 ep=ep+sccdist(i,j)*texp
402 fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
404 C write(*,'(2i8,6f12.5)') i, j, dist,
405 C & fdist(i,j), ep, fp, sccdist(i,j), scc
413 c IF(ABS(EP).lt.1.0d-20)THEN
416 c IF (ABS(FP).lt.1.0d-20) THEN
420 edfadis=edfadis+ep*dis_inc*wwdist
422 gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
423 gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
424 gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
426 gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
427 gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
428 gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
435 subroutine edfat(edfator)
437 implicit real*8 (a-h,o-z)
439 include 'COMMON.CHAIN'
440 include 'COMMON.DERIV'
445 double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
446 double precision cwidth, cwidth2
447 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
448 parameter (TENM20=1.0d-20)
456 do i=idfaphi_start,idfaphi_end
460 iatom(iii)=iphilis(iii,i)+ishiftca
463 C ANGLE VECTOR CALCULTION
464 RIX=C(1,IATOM(2))-C(1,IATOM(1))
465 RIY=C(2,IATOM(2))-C(2,IATOM(1))
466 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
468 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
469 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
470 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
472 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
473 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
474 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
476 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
477 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
478 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
480 GIX=RIY*RIPZ-RIZ*RIPY
481 GIY=RIZ*RIPX-RIX*RIPZ
482 GIZ=RIX*RIPY-RIY*RIPX
484 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
485 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
486 GIPZ=RIPX*RIPPY-RIPY*RIPPX
488 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
489 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
490 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
492 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
493 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
494 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
496 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
497 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
498 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
500 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
501 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
502 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
503 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
505 C END OF ANGLE VECTOR CALCULTION
507 TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
508 APHI(1)=TDOT/(DGI*DRIPP)
509 TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
510 APHI(2)=TDOT/(DGIP*DRIP3)
518 DDPS1=APHI(1)-FPHI1(i,j)
519 DDPS2=APHI(2)-FPHI2(i,j)
521 DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
523 if (dtmp.ge.15.0d0) then
526 c ps_tmp = dfaexp(idint(dtmp*1000)+1)
530 ephi=ephi+sccphi(i,j)*ps_tmp
532 tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
533 tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
536 C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
537 C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
540 ephi=-ephi/scc*phi_inc*wwangle
541 tfphi1=tfphi1/scc*phi_inc*wwangle
542 tfphi2=tfphi2/scc*phi_inc*wwangle
544 IF (ABS(EPHI).LT.1d-20) THEN
547 IF (ABS(TFPHI1).LT.1d-20) THEN
550 IF (ABS(TFPHI2).LT.1d-20) THEN
554 C FORCE DIRECTION CALCULATION
559 DM1=1.0d0/(DGI*DRIPP)
561 GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
562 DM2=GIRPP/(DGI**3*DRIPP)
563 DM3=GIRPP/(DGI*DRIPP**3)
565 DM4=1.0d0/(DGIP*DRIP3)
567 GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
568 DM5=GIRP3/(DGIP**3*DRIP3)
569 DM6=GIRP3/(DGIP*DRIP3**3)
571 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
572 & +( GIZ* RIPY- GIY* RIPZ)*DM2
573 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
574 & +( GIX* RIPZ- GIZ* RIPX)*DM2
575 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
576 & +( GIY* RIPX- GIX* RIPY)*DM2
580 C SECOND ATOM BY PHI1
581 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
582 & -(CIPY*GIZ-CIPZ*GIY)*DM2
583 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
584 & -(CIPZ*GIX-CIPX*GIZ)*DM2
585 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
586 & -(CIPX*GIY-CIPY*GIX)*DM2
590 C SECOND ATOM BY PHI2
592 & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
593 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
595 & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
596 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
598 & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
599 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
601 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
602 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
603 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
604 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
605 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
606 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
612 & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
613 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
615 & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
616 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
618 & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
619 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
620 C FOURTH ATOM BY PHI1
621 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
622 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
623 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
624 C FOURTH ATOM BY PHI2
626 & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
627 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
628 & + RIP3X*DM6)*TFPHI2
630 & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
631 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
632 & + RIP3Y*DM6)*TFPHI2
634 & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
635 & -( GIPX*RIPY-RIPX*GIPY)*DM5
636 & + RIP3Z*DM6)*TFPHI2
638 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
639 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
640 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
641 C END OF FORCE DIRECTION
644 gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
645 gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
646 gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
649 enephi = enephi + ephi
650 c end of single assignment statement
652 C END OF PHI RESTRAINT
654 C START OF THETA ANGLE
655 do i=idfathe_start,idfathe_end
659 iatom(iii)=ithelis(iii,i)+ishiftca
663 C ANGLE VECTOR CALCULTION
664 RIX=C(1,IATOM(2))-C(1,IATOM(1))
665 RIY=C(2,IATOM(2))-C(2,IATOM(1))
666 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
668 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
669 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
670 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
672 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
673 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
674 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
676 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
677 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
678 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
680 GIX=RIY*RIPZ-RIZ*RIPY
681 GIY=RIZ*RIPX-RIX*RIPZ
682 GIZ=RIX*RIPY-RIY*RIPX
684 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
685 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
686 GIPZ=RIPX*RIPPY-RIPY*RIPPX
688 GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
689 GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
690 GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
692 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
693 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
694 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
696 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
697 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
698 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
700 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
701 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
702 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
704 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
705 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
706 DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
707 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
708 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
709 C END OF ANGLE VECTOR CALCULTION
711 TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
712 ATHE(1)=TDOT/(DGI*DGIP)
713 TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
714 ATHE(2)=TDOT/(DGIP*DGIPP)
723 ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
724 ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
725 dtmp= (ddth1**2+ddth2**2)/cwidth2
726 if ( dtmp .ge. 15.0d0) then
729 c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
733 ethe=ethe+sccthe(i,j)*th_tmp
735 tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
736 tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
738 c write(2,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
739 c & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
742 ethe=-ethe/scc*the_inc*wwangle
743 tfthe1=tfthe1/scc*the_inc*wwangle
744 tfthe2=tfthe2/scc*the_inc*wwangle
745 IF (ABS(ETHE).LT.TENM20) THEN
748 IF (ABS(TFTHE1).LT.TENM20) THEN
751 IF (ABS(TFTHE2).LT.TENM20) THEN
760 DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
761 DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
763 DM4=1.0d0/(DGIP*DGIPP)
764 DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
765 DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
767 C FIRST ATOM BY THETA1
768 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
769 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
770 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
771 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
772 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
773 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
774 C SECOND ATOM BY THETA1
775 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
776 & -(CIPY*GIZ-CIPZ*GIY)*DM2
777 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
778 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
779 & -(CIPZ*GIX-CIPX*GIZ)*DM2
780 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
781 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
782 & -(CIPX*GIY-CIPY*GIX)*DM2
783 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
784 C SECOND ATOM BY THETA2
786 & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
787 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
789 & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
790 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
792 & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
793 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
794 C THIRD ATOM BY THETA1
795 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
796 & -(GIY*RIZ-GIZ*RIY)*DM2
797 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
798 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
799 & -(GIZ*RIX-GIX*RIZ)*DM2
800 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
801 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
802 & -(GIX*RIY-GIY*RIX)*DM2
803 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
804 C THIRD ATOM BY THETA2
806 & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
807 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
808 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
810 & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
811 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
812 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
814 & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
815 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
816 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
817 C FOURTH ATOM BY THETA1
818 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
819 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
820 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
821 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
822 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
823 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
824 C FOURTH ATOM BY THETA2
826 & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
827 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
828 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
830 & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
831 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
832 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
834 & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
835 & -(GIPX*RIPY-GIPY*RIPX)*DM5
836 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
837 C FIFTH ATOM BY THETA2
838 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
839 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
840 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
841 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
842 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
843 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
844 C !! END OF FORCE DIRECTION!!!!
846 gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
847 gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
848 gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
851 enethe = enethe + ethe
854 edfator = enephi + enethe
859 subroutine edfan(edfanei)
860 C DFA neighboring CA restraint
861 implicit real*8 (a-h,o-z)
863 include 'COMMON.CHAIN'
864 include 'COMMON.DERIV'
868 integer kshnum, n1atom
870 double precision enenei,tmp_n
871 double precision pai,hpai
872 double precision jix,jiy,jiz,ndiff,snorm_nei
873 double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
874 double precision dr,dr2,half,ntmp,dtmp
876 parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
877 parameter(pai=3.14159265358979323846D0)
878 parameter(hpai=1.5707963267948966D0)
879 parameter(snorm_nei=0.886226925452758D0)
885 c print*, 's1:', s1(:)
886 c print*, 's2:', s2(:)
888 do i=idfanei_start,idfanei_end
891 n1atom=ineilis(i)+ishiftca
892 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
905 do j = ishiftca+1, ilastca
907 if (n1atom.eq.j) cycle
909 jix=c(1,j)-c(1,n1atom)
910 jiy=c(2,j)-c(2,n1atom)
911 jiz=c(3,j)-c(3,n1atom)
912 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
914 c write(*,*) n1atom, j, dist
917 if (dist.lt.s1(kshnum).and.
918 & dist.gt.s2(kshnum-1)) then
922 c write(*,*) 'case1:',tmp_n
931 elseif(dist.ge.s1(kshnum).and.
932 & dist.le.s2(kshnum)) then
934 dnei=(dist-s1(kshnum))/dr2*pai
935 tmp_n=tmp_n + half*(1+cos(dnei))
936 c write(*,*) 'case2:',tmp_n
937 ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
947 elseif(dist.ge.s1(kshnum-1).and.
948 & dist.le.s2(kshnum-1)) then
949 dnei=(dist-s1(kshnum-1))/dr2*pai
950 tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
951 c write(*,*) 'case3:',tmp_n
952 ftmp = hpai*sin(dnei)/dr2/dist
964 elseif(kshnum.eq.1) then
966 if(dist.lt.s1(kshnum))then
969 c write(*,*) 'case4:',tmp_n
977 elseif(dist.ge.s1(kshnum).and.
978 & dist.le.s2(kshnum))then
980 dnei=(dist-s1(kshnum))/dr2*pai
981 tmp_n=tmp_n + half*(1+cos(dnei))
982 c write(*,*) 'case5:',tmp_n
983 ftmp = -hpai*sin(dnei)/dr2/dist
1002 do imin=1,ineinum(i)
1004 ndiff = tmp_n-fnei(i,imin)
1007 if (dtmp.ge.15.0d0) then
1010 c ntmp = dfaexp( idint(dtmp*1000) + 1 )
1014 enei=enei+sccnei(i,imin)*ntmp
1016 & sccnei(i,imin)*ntmp*ndiff*2.0d0
1017 scc=scc+sccnei(i,imin)
1019 c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
1020 c & fnei(i,imin),sccnei(i,imin),enei,scc
1023 enei=-enei/scc*snorm_nei*nei_inc*wwnei
1024 tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
1026 c if (abs(enei).lt.1.0d-20)then
1029 c if (abs(tmp_fnei).lt.1.0d-20) then
1038 do j=ishiftca+1,ilastca
1039 t2dx(j)=t2dx(j)*tmp_fnei
1040 t2dy(j)=t2dy(j)*tmp_fnei
1041 t2dz(j)=t2dz(j)*tmp_fnei
1044 gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
1045 gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
1046 gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
1048 do j=ishiftca+1,ilastca
1049 gdfan(1,j)=gdfan(1,j)+t2dx(j)
1050 gdfan(2,j)=gdfan(2,j)+t2dy(j)
1051 gdfan(3,j)=gdfan(3,j)+t2dz(j)
1053 c energy calculation
1064 subroutine edfab(edfabeta)
1066 implicit real*8 (a-h,o-z)
1068 include 'DIMENSIONS'
1069 include 'COMMON.CHAIN'
1070 include 'COMMON.DERIV'
1071 include 'COMMON.DFA'
1074 parameter(PAI=3.14159265358979323846D0)
1075 parameter (maxca=800)
1077 real*8 bx(maxres),by(maxres),bz(maxres)
1078 real*8 vbet(maxres,maxres)
1079 real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
1080 real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
1081 real*8 vbeta,vbetp,vbetm
1082 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
1083 & c00,s00,ulnex,dnex
1084 real*8 dp45,dm45,w_beta
1086 real*8 cph(maxca),cth(maxca)
1087 real*8 atx(maxca),aty(maxca),atz(maxca)
1088 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1089 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1090 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1092 real*8 astx(maxca),asty(maxca),astz(maxca)
1093 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1094 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1095 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1097 real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
1098 & astxnum(maxca),astynum(maxca),astznum(maxca),
1099 & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
1100 & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
1101 & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
1102 & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
1103 & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
1104 & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
1105 & cth_orig(maxca),sth_orig(maxca)
1107 common /sheca/ bx,by,bz
1108 common /shee/ vbeta,vbet,vbetp,vbetm
1109 common /shetf/ shetfx,shetfy,shetfz
1110 common /shef/ shefx, shefy, shefz
1111 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
1112 & c00,s00,ulnex,dnex
1113 common /sheconst/ dp45,dm45,w_beta
1115 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1116 $ atmmz,atm3x,atm3y,atm3z
1117 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1118 $ astmmz,astm3x,astm3y,astm3z
1120 common /coscos/ cph,cth
1123 C End of sheet variables
1126 double precision enebet
1129 c bx=0.0d0;by=0.0d0;bz=0.0d0
1130 c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1134 do i=ishiftca+1,ilastca
1135 bx(i-ishiftca)=c(1,i)
1136 by(i-ishiftca)=c(2,i)
1137 bz(i-ishiftca)=c(3,i)
1140 c do i=1,ilastca-ishiftca
1141 c read(99,*) bx(i),by(i),bz(i)
1149 ULNEX=COS(60.0D0/180.0D0*PAI)
1156 C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1157 S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1163 C END OF INITIALIZATION
1165 nca=ilastca-ishiftca
1167 call angvectors(nca)
1168 call sheetforce(nca,wshet)
1170 c end of sheet energy and force
1173 shetfx(j)=shetfx(j)*beta_inc
1174 shetfy(j)=shetfy(j)*beta_inc
1175 shetfz(j)=shetfz(j)*beta_inc
1176 c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1179 vbeta=vbeta*beta_inc
1184 gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1185 gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1186 gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1191 write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1205 call angvectors(nca)
1207 call angvectors(nca)
1208 atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1209 astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1211 atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1212 astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1215 atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1216 astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1219 atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1220 astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1224 call angvectors(nca)
1226 call angvectors(nca)
1228 atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1229 astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1231 atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1232 astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1235 atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1236 astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1239 atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1240 astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1244 call angvectors(nca)
1246 call angvectors(nca)
1249 atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1250 astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1252 atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1253 astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1256 atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1257 astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1260 atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1261 astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1267 write (*,'(2i5,a2,6f10.5)')
1268 & i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1269 & astxnum(i),astx(i),astxnum(i)/astx(i),
1270 & i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1271 & astynum(i),asty(i),astynum(i)/asty(i),
1272 & i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1273 & astznum(i),astz(i),astznum(i)/astz(i),
1274 & i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1275 & astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1276 & i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1277 & astmynum(i),astmy(i),astmynum(i)/astmy(i),
1278 & i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1279 & astmznum(i),astmz(i),astmznum(i)/astmz(i),
1280 & i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1281 & astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1282 & i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1283 & astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1284 & i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1285 & astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1286 & i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1287 & astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1288 & i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1289 & astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1290 & i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1291 & astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1292 & i,0," ",cth_orig(i),sth_orig(i)
1302 call angvectors(nca)
1303 call sheetforce(nca,wshet)
1304 vbeta1=vbeta*beta_inc
1306 call angvectors(nca)
1307 call sheetforce(nca,wshet)
1308 vbeta2=vbeta*beta_inc
1309 gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1313 call angvectors(nca)
1314 call sheetforce(nca,wshet)
1315 vbeta1=vbeta*beta_inc
1317 call angvectors(nca)
1318 call sheetforce(nca,wshet)
1319 vbeta2=vbeta*beta_inc
1320 gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1324 call angvectors(nca)
1325 call sheetforce(nca,wshet)
1326 vbeta1=vbeta*beta_inc
1328 call angvectors(nca)
1329 call sheetforce(nca,wshet)
1330 vbeta2=vbeta*beta_inc
1331 gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1338 call angvectors(nca)
1339 call sheetforce(nca,wshet)
1341 shetfx(j)=shetfx(j)*beta_inc
1342 shetfy(j)=shetfy(j)*beta_inc
1343 shetfz(j)=shetfz(j)*beta_inc
1347 write(*,*) 'xyz analytical and numerical gradient'
1349 write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1350 & ,(-gdfab(i,j),i=1,3)
1354 write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1355 & shetfy(j)/gdfab(2,j),
1356 & shetfz(j)/gdfab(3,j)
1364 C-------------------------------------------------------------------------------
1365 subroutine angvectors(nca)
1366 c implicit real*4(a-h,o-z)
1370 parameter(maxca=800)
1372 parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1374 real*8 bx(maxca),by(maxca),bz(maxca)
1375 real*8 dis(maxca,maxca)
1376 real*8 apx(maxca),apy(maxca),apz(maxca)
1377 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
1378 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
1379 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
1380 real*8 atx(maxca),aty(maxca),atz(maxca)
1381 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1382 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1383 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1384 real*8 astx(maxca),asty(maxca),astz(maxca)
1385 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1386 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1387 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1389 real*8 cph(maxca),cth(maxca)
1392 integer i, ip, ipp, ip3, j
1393 real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1394 real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1395 real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1396 real*8 cix, ciy, ciz, cipx, cipy, cipz
1397 real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1398 real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
1399 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1400 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1401 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1402 real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1403 real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1404 real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1405 real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1406 real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1407 real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1408 real*8 grpp,gx,gy,gz
1409 real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1410 real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1411 integer inb,nmax,iselect
1413 common /sheca/ bx,by,bz
1414 common /difvec/ rx, ry, rz
1415 common /ulang/ ulcos
1416 common /phys1/ inb,nmax,iselect
1419 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1420 & apmmz,apm3x,apm3y,apm3z
1421 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1422 & atmmz,atm3x,atm3y,atm3z
1423 common /coscos/ cph,cth
1424 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1425 & astmmz,astm3x,astm3y,astm3z
1427 C-------------------------------------------------------------------------------
1428 c write(*,*) 'inside angvectors'
1433 cph=zero; cth=zero; sth=zero
1434 apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1435 apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1436 atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1437 atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1438 astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1439 astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1442 C r[x,y,z] calc and distance calculation
1443 rx=zero;ry=zero;rz=zero
1450 dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1451 c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1452 c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1453 c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1454 c write(*,*) 'dis(i,j):',i,j,dis(i,j)
1457 c end of r[x,y,z] calc
1463 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1464 ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1465 $ +rz(i,ip)*rz(ip,ipp)
1466 ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1469 c end of virtual bond angle
1470 c write(*,*) 'inside angvectors1'
1482 rippx=bx(ip3)-bx(ipp)
1483 rippy=by(ip3)-by(ipp)
1484 rippz=bz(ip3)-bz(ipp)
1486 gx=riy*ripz-riz*ripy
1487 gy=riz*ripx-rix*ripz
1488 gz=rix*ripy-riy*ripx
1489 gpx=ripy*rippz-ripz*rippy
1490 gpy=ripz*rippx-ripx*rippz
1491 gpz=ripx*rippy-ripy*rippx
1492 gpcrp_x=gpy*ripz-gpz*ripy
1493 gpcrp_y=gpz*ripx-gpx*ripz
1494 gpcrp_z=gpx*ripy-gpy*ripx
1495 d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1496 gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1497 & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1503 gmx=rimy*riz-rimz*riy
1504 gmy=rimz*rix-rimx*riz
1505 gmz=rimx*riy-rimy*rix
1506 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1508 ggm=gmx*gx+gmy*gy+gmz*gz
1509 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1515 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1517 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1518 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1520 c write(*,*) 'inside angvectors2'
1522 rimmx=bx(i-1)-bx(i-2)
1523 rimmy=by(i-1)-by(i-2)
1524 rimmz=bz(i-1)-bz(i-2)
1526 gmmx=rimmy*rimz-rimmz*rimy
1527 gmmy=rimmz*rimx-rimmx*rimz
1528 gmmz=rimmx*rimy-rimmy*rimx
1529 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1531 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1532 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1533 gmcrim_x=gmy*rimz-gmz*rimy
1534 gmcrim_y=gmz*rimx-gmx*rimz
1535 gmcrim_z=gmx*rimy-gmy*rimx
1536 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1537 d_gmcrim3=d_gmcrim**3
1538 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1539 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1543 rim3x=bx(i-2)-bx(i-3)
1544 rim3y=by(i-2)-by(i-3)
1545 rim3z=bz(i-2)-bz(i-3)
1546 g3x=rim3y*rimmz-rim3z*rimmy
1547 g3y=rim3z*rimmx-rim3x*rimmz
1548 g3z=rim3x*rimmy-rim3y*rimmx
1549 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1550 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1551 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1552 cc**********************************************************************
1553 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1554 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1555 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1556 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1557 d_gmmcrimm3=d_gmmcrimm**3
1558 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1559 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1566 dg=sqrt(gx**2+gy**2+gz**2)
1567 dgp=sqrt(gpx**2+gpy**2+gpz**2)
1570 ggp=gx*gpx+gy*gpy+gz*gpz
1571 grpp=gx*rippx+gy*rippy+gz*rippz
1573 if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1574 & .and.d_gpcrp.gt.0.0D0) then
1575 cph(i)=grpp/dg/dripp
1577 sth(i)=gpcrp__g/d_gpcrp/dg
1585 c write(*,*) 'inside angvectors3'
1587 if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1588 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1591 d12=1.0D0/(dg*dripp)
1592 d13=grpp/(dg3*dripp)
1593 sd10=1.0D0/(d_gpcrp*dg)
1594 sd11=gpcrp__g/(d_gpcrp*dg3)
1604 atx(i)=(ripz*gpy-ripy*gpz)*d10
1605 & -(gy*ripz-gz*ripy)*d11
1606 aty(i)=(ripx*gpz-ripz*gpx)*d10
1607 & -(gz*ripx-gx*ripz)*d11
1608 atz(i)=(ripy*gpx-ripx*gpy)*d10
1609 & -(gx*ripy-gy*ripx)*d11
1610 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1611 & +ripy*gpy*ripx-gpx*ripz**2)
1612 & -sd11*(gy*ripz-gz*ripy)
1613 asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1614 & -gpy*ripx**2+gpz*ripy*ripz)
1615 & -sd11*(-gx*ripz+gz*ripx)
1616 astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1617 & -gpz*ripy**2+ripz*gpx*ripx)
1618 & -sd11*(gx*ripy-gy*ripx)
1619 apx(i)=(ripz*rippy-ripy*rippz)*d12
1620 & -(gy*ripz-gz*ripy)*d13
1621 apy(i)=(ripx*rippz-ripz*rippx)*d12
1622 & -(gz*ripx-gx*ripz)*d13
1623 apz(i)=(ripy*rippx-ripx*rippy)*d12
1624 & -(gx*ripy-gy*ripx)*d13
1636 if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1637 & .and.d_gcr3.gt.0.0D0) then
1641 d23=1.0D0/(dgm*drip)
1642 d24=gmrp/(dgm3*drip)
1643 sd20=1.0D0/(d_gcr*dgm)
1644 sd21=gcr__gm/(d_gcr3*dgm)
1645 sd22=gcr__gm/(d_gcr*dgm3)
1656 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1657 & -(ciy*gmz-ciz*gmy)*d21
1658 & +(ripy*gz-ripz*gy)*d22
1659 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1660 & -(ciz*gmx-cix*gmz)*d21
1661 & +(ripz*gx-ripx*gz)*d22
1662 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1663 & -(cix*gmy-ciy*gmx)*d21
1664 & +(ripx*gy-ripy*gx)*d22
1665 cc**********************************************************************
1666 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1667 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1668 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1669 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1670 & +gcr_z*(-ripz*rix+gy))
1671 & -sd22*(-gmy*ciz+gmz*ciy)
1673 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1674 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1676 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1677 & -gcr_z*(ripz*riy+gx))
1678 & -sd22*(gmx*ciz-gmz*cix)
1680 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1681 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1683 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1684 & +gcr_z*(ripy*riy+ripx*rix))
1685 & -sd22*(-gmx*ciy+gmy*cix)
1686 cc**********************************************************************
1687 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1688 & -(ciy*gmz-ciz*gmy)*d24
1689 apmy(i)=(ciz*ripx-ripz*cix)*d23
1690 & -(ciz*gmx-cix*gmz)*d24
1691 apmz(i)=(cix*ripy-ripx*ciy)*d23
1692 & -(cix*gmy-ciy*gmx)*d24
1696 if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1697 & .and.d_gmcrim3.gt.0.0D0) then
1698 d30=1.0D0/(dgm*dgmm)
1699 d31=gmmgm/(dgm3*dgmm)
1700 d32=gmmgm/(dgm*dgmm3)
1701 d33=1.0D0/(dgmm*dri)
1702 d34=gmmr/(dgmm3*dri)
1703 d35=gmmr/(dgmm*dri3)
1704 sd30=1.0D0/(d_gmcrim*dgmm)
1705 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1706 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1719 c write(*,*) 'inside angvectors4'
1721 cc**********************************************************************
1722 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1723 & -(ciy*gmz-ciz*gmy)*d31
1724 & -(gmmy*rimmz-gmmz*rimmy)*d32
1725 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1726 & -(ciz*gmx-cix*gmz)*d31
1727 & -(gmmz*rimmx-gmmx*rimmz)*d32
1728 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1729 & -(cix*gmy-ciy*gmx)*d31
1730 & -(gmmx*rimmy-gmmy*rimmx)*d32
1731 cc**********************************************************************
1732 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1733 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1734 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1735 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1736 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1737 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1738 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1740 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1741 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1742 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1743 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1744 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1745 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1746 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1748 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1749 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1750 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1751 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1752 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1753 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1754 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1755 c**********************************************************************
1756 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1757 & -(gmmy*rimmz-gmmz*rimmy)*d34
1759 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1760 & -(gmmz*rimmx-gmmx*rimmz)*d34
1762 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1763 & -(gmmx*rimmy-gmmy*rimmx)*d34
1768 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1769 & .and.drim3.gt.0.0D0
1770 & .and.d_gmmcrimm3.gt.0.0D0) then
1771 d40=1.0D0/(dg30*dgmm)
1772 d41=g3gmm/(dg30*dgmm3)
1773 d42=1.0D0/(dg30*drim)
1774 d43=g3rim_/(dg30*drim3)
1775 sd40=1.0D0/(dg30*d_gmmcrimm)
1776 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1785 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1786 & -(gmmy*rimmz-gmmz*rimmy)*d41
1787 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1788 & -(gmmz*rimmx-gmmx*rimmz)*d41
1789 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1790 & -(gmmx*rimmy-gmmy*rimmx)*d41
1791 cc**********************************************************************
1792 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1793 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1794 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1795 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1797 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1798 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1799 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1800 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1802 c & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1804 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1805 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1806 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1807 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1808 c**********************************************************************
1809 apm3x(i)=g3x*d42-rimx*d43
1810 apm3y(i)=g3y*d42-rimy*d43
1811 apm3z(i)=g3z*d42-rimz*d43
1814 c*******************************************************************************
1816 c write(*,*) 'inside angvectors5'
1823 rimmx=bx(i-1)-bx(i-2)
1824 rimmy=by(i-1)-by(i-2)
1825 rimmz=bz(i-1)-bz(i-2)
1826 rim3x=bx(i-2)-bx(i-3)
1827 rim3y=by(i-2)-by(i-3)
1828 rim3z=bz(i-2)-bz(i-3)
1829 gmmx=rimmy*rimz-rimmz*rimy
1830 gmmy=rimmz*rimx-rimmx*rimz
1831 gmmz=rimmx*rimy-rimmy*rimx
1832 g3x=rim3y*rimmz-rim3z*rimmy
1833 g3y=rim3z*rimmx-rim3x*rimmz
1834 g3z=rim3x*rimmy-rim3y*rimmx
1836 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1837 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1838 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1843 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1844 cc**********************************************************************
1845 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1846 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1847 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1848 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1849 d_gmmcrimm3=d_gmmcrimm**3
1850 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1851 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1853 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1854 & .and.drim3.gt.0.0D0
1855 & .and.d_gmmcrimm3.gt.0.0D0) then
1856 d40=1.0D0/(dg30*dgmm)
1857 d41=g3gmm/(dg30*dgmm3)
1858 d42=1.0D0/(dg30*drim)
1859 d43=g3rim_/(dg30*drim3)
1860 sd40=1.0D0/(dg30*d_gmmcrimm)
1861 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1870 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1871 & -(gmmy*rimmz-gmmz*rimmy)*d41
1872 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1873 & -(gmmz*rimmx-gmmx*rimmz)*d41
1874 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1875 & -(gmmx*rimmy-gmmy*rimmx)*d41
1876 cc**********************************************************************
1877 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1878 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1879 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1880 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1882 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1883 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1884 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1885 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1887 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1888 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1889 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1890 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1891 cc**********************************************************************
1892 apm3x(i)=g3x*d42-rimx*d43
1893 apm3y(i)=g3y*d42-rimy*d43
1894 apm3z(i)=g3z*d42-rimz*d43
1904 gmx=rimy*riz-rimz*riy
1905 gmy=rimz*rix-rimx*riz
1906 gmz=rimx*riy-rimy*rix
1907 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1911 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1912 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1913 gmcrim_x=gmy*rimz-gmz*rimy
1914 gmcrim_y=gmz*rimx-gmx*rimz
1915 gmcrim_z=gmx*rimy-gmy*rimx
1916 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1917 d_gmcrim3=d_gmcrim**3
1918 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1919 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1921 if(dgm3.gt.0.0D0.and.
1922 & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1923 & .and.d_gmcrim3.gt.0.0D0) then
1924 d30=1.0D0/(dgm*dgmm)
1925 d31=gmmgm/(dgm3*dgmm)
1926 d32=gmmgm/(dgm*dgmm3)
1927 d33=1.0D0/(dgmm*dri)
1928 d34=gmmr/(dgmm3*dri)
1929 d35=gmmr/(dgmm*dri3)
1930 sd30=1.0D0/(d_gmcrim*dgmm)
1931 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1932 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1945 cc**********************************************************************
1946 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1947 & -(ciy*gmz-ciz*gmy)*d31
1948 & -(gmmy*rimmz-gmmz*rimmy)*d32
1949 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1950 & -(ciz*gmx-cix*gmz)*d31
1951 & -(gmmz*rimmx-gmmx*rimmz)*d32
1952 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1953 & -(cix*gmy-ciy*gmx)*d31
1954 & -(gmmx*rimmy-gmmy*rimmx)*d32
1955 cc**********************************************************************
1956 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1957 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1958 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1959 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1960 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1961 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1962 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1964 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1965 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1966 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1967 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1968 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1969 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1970 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1972 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1973 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1974 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1975 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1976 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1977 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1978 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1979 cc**********************************************************************
1980 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1981 & -(gmmy*rimmz-gmmz*rimmy)*d34
1983 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1984 & -(gmmz*rimmx-gmmx*rimmz)*d34
1986 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1987 & -(gmmx*rimmy-gmmy*rimmx)*d34
1991 c write(*,*) 'inside angvectors6'
2001 gx=riy*ripz-riz*ripy
2002 gy=riz*ripx-rix*ripz
2003 gz=rix*ripy-riy*ripx
2004 ggm=gmx*gx+gmy*gy+gmz*gz
2005 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
2006 dg=sqrt(gx**2+gy**2+gz**2)
2012 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
2014 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
2015 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
2016 if(dgm3.gt.0.0D0.and.
2017 & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
2022 d23=1.0D0/(dgm*drip)
2023 d24=gmrp/(dgm3*drip)
2024 sd20=1.0D0/(d_gcr*dgm)
2025 sd21=gcr__gm/(d_gcr3*dgm)
2026 sd22=gcr__gm/(d_gcr*dgm3)
2037 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
2038 & -(ciy*gmz-ciz*gmy)*d21
2039 & +(ripy*gz-ripz*gy)*d22
2040 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
2041 & -(ciz*gmx-cix*gmz)*d21
2042 & +(ripz*gx-ripx*gz)*d22
2043 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
2044 & -(cix*gmy-ciy*gmx)*d21
2045 & +(ripx*gy-ripy*gx)*d22
2046 cc**********************************************************************
2047 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
2048 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
2049 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
2050 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
2051 & +gcr_z*(-ripz*rix+gy))
2052 & -sd22*(-gmy*ciz+gmz*ciy)
2054 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
2055 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
2057 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
2058 & -gcr_z*(ripz*riy+gx))
2059 & -sd22*(gmx*ciz-gmz*cix)
2061 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
2062 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
2064 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
2065 & +gcr_z*(ripy*riy+ripx*rix))
2066 & -sd22*(-gmx*ciy+gmy*cix)
2067 cc**********************************************************************
2069 apmx(i)=(ciy*ripz-ripy*ciz)*d23
2070 & -(ciy*gmz-ciz*gmy)*d24
2071 apmy(i)=(ciz*ripx-ripz*cix)*d23
2072 & -(ciz*gmx-cix*gmz)*d24
2073 apmz(i)=(cix*ripy-ripx*ciy)*d23
2074 & -(cix*gmy-ciy*gmx)*d24
2082 c-------------------------------------------------------------------------------
2083 C---------------------------------------------------------------------------------
2084 subroutine sheetforce(nca,wshet)
2087 c this should be matched with dfa.fcm
2089 parameter(maxca=800)
2090 cc**********************************************************************
2093 integer inb,nmax,iselect
2095 c real*8 dfaexp(15001)
2097 real*8 vbeta,vbetp,vbetm
2098 real*8 shefx(maxca,12)
2099 real*8 shefy(maxca,12),shefz(maxca,12)
2100 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2101 real*8 vbet(maxca,maxca)
2102 real*8 wshet(maxca,maxca)
2103 real*8 bx(maxca),by(maxca),bz(maxca)
2105 common /sheca/ bx,by,bz
2106 common /phys1/ inb,nmax,iselect
2107 common /shef/ shefx,shefy,shefz
2108 common /shee/ vbeta,vbet,vbetp,vbetm
2109 common /shetf/ shetfx,shetfy,shetfz
2126 call sheetene(nca,wshet)
2129 887 format(a,1x,i6,3x,f12.8)
2130 888 format(a,1x,i4,1x,i4,3x,f12.8)
2131 889 format(a,1x,i4,3x,f12.8)
2132 !write(2,*) 'coord : '
2134 !write(2,887) 'bx:',i,bx(i)
2135 !write(2,887) 'by:',i,by(i)
2136 !write(2,887) 'bz:',i,bz(i)
2138 !write(2,*) 'After sheetforce1'
2141 !write(2,888) 'shefx :',i,k,shefx(i,k)
2142 !write(2,888) 'shefy :',i,k,shefy(i,k)
2143 !write(2,888) 'shefz :',i,k,shefz(i,k)
2149 !write(2,*) 'After sheetforce5'
2152 !write(2,888) 'shefx :',i,k,shefx(i,k)
2153 !write(2,888) 'shefy :',i,k,shefy(i,k)
2154 !write(2,888) 'shefz :',i,k,shefz(i,k)
2160 !write(2,*) 'After sheetforce6'
2163 !write(2,888) 'shefx :',i,k,shefx(i,k)
2164 !write(2,888) 'shefy :',i,k,shefy(i,k)
2165 !write(2,888) 'shefz :',i,k,shefz(i,k)
2171 !write(2,*) 'After sheetforce11'
2174 !write(2,888) 'shefx :',i,k,shefx(i,k)
2175 !write(2,888) 'shefy :',i,k,shefy(i,k)
2176 !write(2,888) 'shefz :',i,k,shefz(i,k)
2182 !write(2,*) 'After sheetforce12'
2185 !write(2,888) 'shefx :',i,k,shefx(i,k)
2186 !write(2,888) 'shefy :',i,k,shefy(i,k)
2187 !write(2,888) 'shefz :',i,k,shefz(i,k)
2193 shetfx(i)=shetfx(i)+shefx(i,k)
2194 shetfy(i)=shetfy(i)+shefy(i,k)
2195 shetfz(i)=shetfz(i)+shefz(i,k)
2198 !write(2,*) 'Beta Finished'
2200 !write(2,889) 'shetfx : ',i,shetfx(i)
2201 !write(2,889) 'shetfy : ',i,shetfy(i)
2202 !write(2,889) 'shetfz : ',i,shetfz(i)
2208 c-------------------------------------------------------------------------------
2209 subroutine sheetene(nca,wshet)
2212 parameter(maxca=800)
2213 real*8 dfa_cutoff,dfa_cutoff_delta
2214 parameter(dfa_cutoff=15.5d0)
2215 parameter(dfa_cutoff_delta=0.5d0)
2216 cc******************************************************************************
2218 c real*8 dfaexp(15001)
2219 real*8 dtmp1, dtmp2, dtmp3
2221 real*8 vbet(maxca,maxca)
2222 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2223 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2224 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2225 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2226 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2227 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2228 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2229 real*8 cph(maxca),cth(maxca)
2230 real*8 rx(maxca,maxca)
2231 real*8 ry(maxca,maxca),rz(maxca,maxca)
2232 real*8 bx(maxca),by(maxca),bz(maxca)
2233 real*8 dis(maxca,maxca)
2235 cc**********************************************************************
2236 real*8 astx(maxca),asty(maxca),astz(maxca)
2237 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2238 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2239 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2241 real*8 wshet(maxca,maxca)
2242 real*8 dp45, dm45, w_beta
2243 real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2245 integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2247 real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2249 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2250 common /shetf/ shetfx,shetfy,shetfz
2252 common /sheca/ bx,by,bz
2253 common /phys1/ inb,nmax,iselect
2255 common /difvec/ rx,ry,rz
2256 common /coscos/ cph,cth
2257 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2258 & c00,s00,ulnex,dnex
2259 common /sheconst/ dp45,dm45,w_beta
2260 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2261 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2262 common /shee/ vbeta,vbet,vbetp,vbetm
2263 common /ulang/ ulcos
2264 cc**********************************************************************
2265 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2266 & astmmz,astm3x,astm3y,astm3z
2269 real*8 r_pair_mat(maxca,maxca)
2270 real*8 e_gcont,fprim_gcont,de_gcont
2271 ci integer istrand(maxca,maxca)
2272 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2273 ci common /shetest/ istrand,istrand_p,istrand_m
2274 common /beta_p/ r_pair_mat
2275 C-------------------------------------------------------------------------------
2279 r_pair_mat(i,j)=wshet(i,j)
2280 c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2292 if (dis(i,j).lt.dfa_cutoff) then
2293 call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2294 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2301 cc**********************************************************************
2302 y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2303 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2305 y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2306 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2319 ci if(istrand(i,j).eq.1) then
2321 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2322 yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2325 pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2326 $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2327 pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2328 $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2329 pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2330 $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2331 pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2332 $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2334 yshe1=pin1(i,j)**2+pin2(i,j)**2
2335 yshe1=-0.5d0*yshe1/dshe
2336 yshe2=pin3(i,j)**2+pin4(i,j)**2
2337 yshe2=-0.5d0*yshe2/dshe
2339 ci if((yshe1+yshe2).ge.-4) then
2346 C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2347 C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2348 C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2349 C write(*,*) 'dis(i,j):',i,j,dis(i,j)
2350 C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2351 C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2352 C write(*,*) 'pin1:',pin1(i,j)
2353 C write(*,*) 'pin2:',pin2(i,j)
2354 C write(*,*) 'pin3:',pin3(i,j)
2355 C write(*,*) 'pin4:',pin4(i,j)
2358 C write(*,*) 'yy1:',yy1
2359 C write(*,*) 'yy2:',yy2
2360 C write(*,*) 'yshe1:',yshe1
2361 C write(*,*) 'yshe2:',yshe2
2364 ci if (istrand_p(i,j).eq.1) then
2372 dtmp3 = y+yy1+yy2+yshe1+yshe2
2374 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2375 C write(*,*)'2', y,yy1,yy2
2376 C write(*,*)'3', yshe1,yshe2
2378 cc if (dtmp3.le.-35.0d0) then
2379 c vbetap(i,j)=-dp45*exp(dtmp3)
2380 cc vbetap(i,j)=0.0d0
2382 c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2383 vbetap(i,j)=-dp45*exp(dtmp3)
2386 cc if (dtmp1.le.-35.0d0) then
2387 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2388 cc vbetap1(i,j)=0.0d0
2390 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2391 c $ *dfaexp(idint(-dtmp1*1000)+1)
2392 vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2395 cc if (dtmp2.le.-35.0d0) then
2396 C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2397 cc vbetap2(i,j)=0.0d0
2399 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2400 c $ *dfaexp(idint(-dtmp2*1000)+1)
2401 vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2404 c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2405 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2406 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2408 ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2409 ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2411 ci elseif (istrand_p(i,j).eq.0)then
2417 yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2418 yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2420 pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2421 $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2422 pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2423 $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2424 pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2425 $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2426 pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2427 $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2429 yshe1=pina1(i,j)**2+pina2(i,j)**2
2430 yshe1=-0.5d0*yshe1/dshe
2431 yshe2=pina3(i,j)**2+pina4(i,j)**2
2432 yshe2=-0.5d0*yshe2/dshe
2434 ci if((yshe1+yshe2).ge.-4) then
2441 C write(*,*) 'pina1:',pina1(i,j)
2442 C write(*,*) 'pina2:',pina2(i,j)
2443 C write(*,*) 'pina3:',pina3(i,j)
2444 C write(*,*) 'pina4:',pina4(i,j)
2445 C write(*,*) 'yshe1:',yshe1
2446 C write(*,*) 'yshe2:',yshe2
2447 C write(*,*) 'dshe:',dshe
2449 ci if (istrand_m(i,j).eq.1) then
2456 dtmp3=y+yy1+yy2+yshe1+yshe2
2460 cc if(dtmp3 .le. -35.0d0) then
2461 c vbetam(i,j)=-dm45*exp(dtmp3)
2462 cc vbetam(i,j)=0.0d0
2464 c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2465 vbetam(i,j)=-dm45*exp(dtmp3)
2468 cc if(dtmp1 .le. -35.0d0) then
2469 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2470 cc vbetam1(i,j)=0.0d0
2472 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2473 c $ *dfaexp(idint(-dtmp1*1000)+1)
2474 vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2477 cc if(dtmp2.le.-35.0d0) then
2478 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2479 cc vbetam2(i,j)=0.0d0
2481 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2482 c $ *dfaexp(idint(-dtmp2*1000)+1)
2483 vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2486 ci elseif (istrand_m(i,j).eq.0)then
2493 c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2494 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2495 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2497 ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2498 ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2500 uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2501 uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2503 c write(*,*) 'uup,uum:', uup, uum
2505 c uup=vbetap1(i,j)+vbetap2(i,j)
2506 c uum=vbetam1(i,j)+vbetam2(i,j)
2511 vbeta=vbeta+vbet(i,j)*e_gcont
2514 if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
2515 c gradient correction from gcont
2516 de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
2517 shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
2518 shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
2519 shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
2521 shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
2522 shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
2523 shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
2525 c energy correction from gcont
2526 vbet(i,j)=vbet(i,j)*e_gcont
2527 vbetap(i,j)=vbetap(i,j)*e_gcont
2528 vbetap1(i,j)=vbetap1(i,j)*e_gcont
2529 vbetap2(i,j)=vbetap2(i,j)*e_gcont
2530 vbetam(i,j)=vbetam(i,j)*e_gcont
2531 vbetam1(i,j)=vbetam1(i,j)*e_gcont
2532 vbetam2(i,j)=vbetam2(i,j)*e_gcont
2536 ci elseif(istrand(i,j).eq.0)then
2540 c write(*,*) 'uup,uum:',uup,uum
2541 c write(*,*) 'vbetap(i,j):',vbetap(i,j)
2542 c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2543 c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2544 c write(*,*) 'vbetam(i,j):',vbetam(i,j)
2545 c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2546 c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2547 c write(*,*) 'uup:',uup
2548 c write(*,*) 'uum:',uum
2549 c write(*,*) 'vbetp:',vbetp
2550 c write(*,*) 'vbetm:',vbetm
2551 c write(*,*) 'vbet(i,j):',vbet(i,j)
2568 ! write(*,*) 'I,J:', i,j
2569 ! write(*,*) 'vbetap(i,j):',vbetap(i,j)
2570 ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2571 ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2572 ! write(*,*) 'vbetam(i,j):',vbetam(i,j)
2573 ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2574 ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2575 ! write(*,*) 'vbet(i,j):',vbet(i,j)
2581 c-------------------------------------------------------------------------------
2582 subroutine sheetforce1
2585 parameter(maxca=800)
2586 real*8 dfa_cutoff,dfa_cutoff_delta
2587 parameter(dfa_cutoff=15.5d0)
2588 parameter(dfa_cutoff_delta=0.5d0)
2589 cc**********************************************************************
2590 real*8 vbet(maxca,maxca)
2591 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2592 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2593 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2594 real*8 cph(maxca),cth(maxca)
2595 real*8 rx(maxca,maxca)
2596 real*8 ry(maxca,maxca),rz(maxca,maxca)
2597 real*8 bx(maxca),by(maxca),bz(maxca)
2598 real*8 dis(maxca,maxca)
2599 real*8 shefx(maxca,12)
2600 real*8 shefy(maxca,12),shefz(maxca,12)
2601 real*8 atx(maxca),aty(maxca),atz(maxca)
2602 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2603 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2604 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2605 real*8 apx(maxca),apy(maxca),apz(maxca)
2606 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2607 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2608 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2610 real*8 astx(maxca),asty(maxca),astz(maxca)
2611 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2612 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2613 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2615 real*8 w_beta,dp45, dm45
2616 real*8 vbeta, vbetp, vbetm
2617 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2618 $ c00,s00,ulnex,dnex
2619 integer inb,nmax,iselect
2621 common /phys1/ inb,nmax,iselect
2623 common /difvec/ rx,ry,rz
2624 common /coscos/ cph,cth
2625 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2626 $ c00,s00,ulnex,dnex
2627 common /sheconst/ dp45,dm45,w_beta
2628 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2629 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2630 $ atmmz,atm3x,atm3y,atm3z
2631 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2632 $ apmmz,apm3x,apm3y,apm3z
2633 common /shef/ shefx,shefy,shefz
2634 common /shee/ vbeta,vbet,vbetp,vbetm
2635 common /ulang/ ulcos
2636 c c**********************************************************************
2637 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2638 $ astmmz,astm3x,astm3y,astm3z
2640 C--------------------------------------------------------------------------------
2642 integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2643 real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2644 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2645 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2646 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2647 real*8 e_gcont,fprim_gcont
2648 C--------------------------------------------------------------------------------
2653 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2658 cc1=(ulcos(imm)-ulnex)/dnex
2659 dmm=cc1/(dis(imm,im)*dis(im,i))
2660 dmm__=cc1*ulcos(imm)/dis(im,i)**2
2661 fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2662 fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2663 fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2667 fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2668 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2669 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2679 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2684 cc1=(ulcos(imm)-ulnex)/dnex
2685 cc2=(ulcos(im)-ulnex)/dnex
2686 dmm1=cc1/(dis(imm,im)*dis(im,i))
2687 dmm2=cc2/(dis(im,i)*dis(i,ip))
2688 dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2689 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2690 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2691 cc**********************************************************************
2692 fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2693 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2694 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2695 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2696 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2697 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2701 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2702 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2703 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2712 c3=(cth(im)*c00+sth(im)*s00-1)/dca
2717 cc2=(ulcos(im)-ulnex)/dnex
2718 cc3=(ulcos(i)-ulnex)/dnex
2719 dmm2=cc2/(dis(im,i)*dis(i,ip))
2720 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2721 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2722 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2723 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2724 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2725 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2726 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2727 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2728 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2729 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2733 fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2734 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2735 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2743 c4=(cth(i)*c00+sth(i)*s00-1)/dca
2748 cc3=(ulcos(i)-ulnex)/dnex
2749 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2750 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2751 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2752 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2753 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2757 fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2758 fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2759 fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2768 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2773 cc7=(ulcos(jmm)-ulnex)/dnex
2774 dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2775 dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2776 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2777 fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2778 fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2782 fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2783 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2784 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2793 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2798 cc7=(ulcos(jmm)-ulnex)/dnex
2799 cc8=(ulcos(jm)-ulnex)/dnex
2800 dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2801 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2802 dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2803 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2804 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2805 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2806 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2807 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2808 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2809 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2810 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2814 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2815 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2816 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2826 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2831 cc8=(ulcos(jm)-ulnex)/dnex
2832 cc9=(ulcos(j)-ulnex)/dnex
2833 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2834 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2835 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2836 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2837 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2838 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2839 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2840 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2841 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2842 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2843 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2847 fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2848 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2849 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2858 c10=(cth(j)*c00+sth(j)*s00-1)/dca
2863 cc9=(ulcos(j)-ulnex)/dnex
2864 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2865 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2866 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2867 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2868 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2872 fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2873 fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2874 fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2882 c----------------------------------------------------------------------------
2883 subroutine sheetforce5
2886 parameter(maxca=800)
2887 real*8 dfa_cutoff,dfa_cutoff_delta
2888 parameter(dfa_cutoff=15.5d0)
2889 parameter(dfa_cutoff_delta=0.5d0)
2890 cc**********************************************************************
2891 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2892 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2893 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2894 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2895 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2896 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2897 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2898 real*8 rx(maxca,maxca)
2899 real*8 ry(maxca,maxca),rz(maxca,maxca)
2900 real*8 bx(maxca),by(maxca),bz(maxca)
2901 real*8 dis(maxca,maxca)
2902 real*8 shefx(maxca,12),shefy(maxca,12)
2903 real*8 shefz(maxca,12)
2904 real*8 dp45,dm45,w_beta
2905 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2906 $ c00,s00,ulnex,dnex
2907 integer inb,nmax,iselect
2908 cc**********************************************************************
2909 common /phys1/ inb,nmax,iselect
2911 common /difvec/ rx,ry,rz
2912 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2913 $ c00,s00,ulnex,dnex
2914 common /sheconst/ dp45,dm45,w_beta
2915 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2916 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2917 common /shef/ shefx,shefy,shefz
2918 ci integer istrand(maxca,maxca)
2919 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2920 ci common /shetest/ istrand,istrand_p,istrand_m
2921 c********************************************************************************
2923 integer i,imm,im,jp,jpp,j
2924 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2925 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2926 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2927 real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2928 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2929 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2930 real*8 e_gcont,fprim_gcont
2931 c********************************************************************************
2937 if (dis(imm,j).lt.dfa_cutoff) then
2938 call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2939 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2944 ci if(istrand(imm,j).eq.1
2945 ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2948 yy1=-(dis(i,jpp)-ulhb)/dlhb
2949 y1x=rx(jpp,i)/dis(i,jpp)
2950 y1y=ry(jpp,i)/dis(i,jpp)
2951 y1z=rz(jpp,i)/dis(i,jpp)
2956 yy33=1.0D0/(dis(im,jp)*dis(im,i))
2957 yyy3=pin1(imm,j)/(dis(im,i)**2)
2958 yy3=-pin1(imm,j)/dshe
2959 y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2960 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2961 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2963 yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2964 yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2965 yyy4b=pin3(imm,j)/(dis(im,i)**2)
2966 yy4=-pin3(imm,j)/dshe
2967 y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2968 $ -yyy4b*rx(im,i))*yy4
2969 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2970 $ -yyy4b*ry(im,i))*yy4
2971 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2972 $ -yyy4b*rz(im,i))*yy4
2975 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2976 yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2977 yy5=-pin4(imm,j)/dshe
2978 y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2979 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2980 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2993 shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2994 $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2995 shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2996 $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2997 shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2998 $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
3000 ! shefx(i,5)=shefx(i,5)
3001 ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
3002 ! shefy(i,5)=shefy(i,5)
3003 ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
3004 ! shefz(i,5)=shefz(i,5)
3005 ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
3007 yy6=-(dis(i,jp)-uldhb)/dldhb
3008 y6x=rx(jp,i)/dis(i,jp)
3009 y6y=ry(jp,i)/dis(i,jp)
3010 y6z=rz(jp,i)/dis(i,jp)
3015 yy88=1.0D0/(dis(im,jpp)*dis(im,i))
3016 yyy8=pina1(imm,j)/(dis(im,i)**2)
3017 yy8=-pina1(imm,j)/dshe
3018 y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
3019 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
3020 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
3022 yy99=1.0D0/(dis(jp,i)*dis(im,i))
3023 yyy9a=pina3(imm,j)/(dis(jp,i)**2)
3024 yyy9b=pina3(imm,j)/(dis(im,i)**2)
3025 yy9=-pina3(imm,j)/dshe
3026 y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
3027 $ -yyy9b*rx(im,i))*yy9
3028 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
3029 $ -yyy9b*ry(im,i))*yy9
3030 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
3031 $ -yyy9b*rz(im,i))*yy9
3033 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
3034 yyy10=pina4(imm,j)/(dis(jp,i)**2)
3035 yy10=-pina4(imm,j)/dshe
3036 y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
3037 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
3038 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
3040 sx=y66x+y8x+y9x+y10x
3041 sy=y66y+y8y+y9y+y10y
3042 sz=y66z+y8z+y9z+y10z
3051 shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
3052 $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3053 shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
3054 $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3055 shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
3056 $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3058 ! shefx(i,5)=shefx(i,5)
3059 ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3060 ! shefy(i,5)=shefy(i,5)
3061 ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3062 ! shefz(i,5)=shefz(i,5)
3063 ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3073 c--------------------------------------------------------------------------c
3074 subroutine sheetforce6
3077 parameter(maxca=800)
3078 real*8 dfa_cutoff,dfa_cutoff_delta
3079 parameter(dfa_cutoff=15.5d0)
3080 parameter(dfa_cutoff_delta=0.5d0)
3081 cc**********************************************************************
3082 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3083 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3084 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3085 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3086 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3087 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3088 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3089 real*8 rx(maxca,maxca)
3090 real*8 ry(maxca,maxca),rz(maxca,maxca)
3091 real*8 bx(maxca),by(maxca),bz(maxca)
3092 real*8 dis(maxca,maxca)
3093 real*8 shefx(maxca,12),shefy(maxca,12)
3094 real*8 shefz(maxca,12)
3095 real*8 dp45,dm45,w_beta
3096 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3097 $ c00,s00,ulnex,dnex
3098 integer inb,nmax,iselect
3099 cc**********************************************************************
3100 common /phys1/ inb,nmax,iselect
3102 common /difvec/ rx,ry,rz
3103 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3104 $ c00,s00,ulnex,dnex
3105 common /sheconst/ dp45,dm45,w_beta
3106 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3107 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3108 common /shef/ shefx,shefy,shefz
3109 ci integer istrand(maxca,maxca)
3110 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3111 ci common /shetest/ istrand,istrand_p,istrand_m
3112 cc**********************************************************************
3114 integer i,imm,im,jp,jpp,j,ip
3115 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3116 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3117 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3118 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3119 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
3120 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
3121 real*8 e_gcont,fprim_gcont
3122 C********************************************************************************
3128 if (dis(im,j).lt.dfa_cutoff) then
3129 call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3130 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3135 ci if(istrand(im,j).eq.1
3136 ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
3139 yy1=-(dis(i,jp)-ulhb)/dlhb
3140 y1x=rx(jp,i)/dis(i,jp)
3141 y1y=ry(jp,i)/dis(i,jp)
3142 y1z=rz(jp,i)/dis(i,jp)
3147 yy33=1.0D0/(dis(i,jp)*dis(i,ip))
3148 yyy3a=pin1(im,j)/(dis(i,jp)**2)
3149 yyy3b=pin1(im,j)/(dis(i,ip)**2)
3150 yy3=-pin1(im,j)/dshe
3151 y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
3152 $ +yyy3b*rx(i,ip))*yy3
3153 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
3154 $ +yyy3b*ry(i,ip))*yy3
3155 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
3156 $ +yyy3b*rz(i,ip))*yy3
3158 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3159 yyy4=pin2(im,j)/(dis(i,jp)**2)
3160 yy4=-pin2(im,j)/dshe
3161 y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3162 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3163 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3165 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3166 yyy5=pin3(im,j)/(dis(i,ip)**2)
3167 yy5=-pin3(im,j)/dshe
3168 y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3169 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3170 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3183 shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3184 $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3185 shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3186 $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3187 shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3188 $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3189 ! shefx(i,6)=shefx(i,6)
3190 ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3191 ! shefy(i,6)=shefy(i,6)
3192 ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3193 ! shefz(i,6)=shefz(i,6)
3194 ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3196 yy6=-(dis(jpp,i)-uldhb)/dldhb
3197 y6x=rx(jpp,i)/dis(jpp,i)
3198 y6y=ry(jpp,i)/dis(jpp,i)
3199 y6z=rz(jpp,i)/dis(jpp,i)
3204 yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3205 yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3206 yyy8b=pina1(im,j)/(dis(i,ip)**2)
3207 yy8=-pina1(im,j)/dshe
3208 y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3209 $ +yyy8b*rx(i,ip))*yy8
3210 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3211 $ +yyy8b*ry(i,ip))*yy8
3212 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3213 $ +yyy8b*rz(i,ip))*yy8
3215 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3216 yyy9=pina2(im,j)/(dis(i,jpp)**2)
3217 yy9=-pina2(im,j)/dshe
3218 y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3219 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3220 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3222 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3223 yyy10=pina3(im,j)/(dis(i,ip)**2)
3224 yy10=-pina3(im,j)/dshe
3225 y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3226 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3227 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3229 sx=y66x+y8x+y9x+y10x
3230 sy=y66y+y8y+y9y+y10y
3231 sz=y66z+y8z+y9z+y10z
3240 shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3241 $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3242 shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3243 $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3244 shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3245 $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3247 ! shefx(i,6)=shefx(i,6)
3248 ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3249 ! shefy(i,6)=shefy(i,6)
3250 ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3251 ! shefz(i,6)=shefz(i,6)
3252 ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3262 c-----------------------------------------------------------------------
3263 subroutine sheetforce11
3266 parameter(maxca=800)
3267 real*8 dfa_cutoff,dfa_cutoff_delta
3268 parameter(dfa_cutoff=15.5d0)
3269 parameter(dfa_cutoff_delta=0.5d0)
3270 cc**********************************************************************
3271 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3272 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3273 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3274 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3275 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3276 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3277 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3278 real*8 rx(maxca,maxca)
3279 real*8 ry(maxca,maxca),rz(maxca,maxca)
3280 real*8 bx(maxca),by(maxca),bz(maxca)
3281 real*8 dis(maxca,maxca)
3282 real*8 shefx(maxca,12),shefy(maxca,12)
3283 real*8 shefz(maxca,12)
3284 real*8 dp45,dm45,w_beta
3285 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3286 $ c00,s00,ulnex,dnex
3287 integer inb,nmax,iselect
3288 cc**********************************************************************
3289 common /phys1/ inb,nmax,iselect
3291 common /difvec/ rx,ry,rz
3292 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3293 $ c00,s00,ulnex,dnex
3294 common /sheconst/ dp45,dm45,w_beta
3295 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3296 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3297 common /shef/ shefx,shefy,shefz
3298 ci integer istrand(maxca,maxca)
3299 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3300 ci common /shetest/ istrand,istrand_p,istrand_m
3301 C********************************************************************************
3303 integer j,jm,jmm,ip,i,ipp
3304 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3305 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3306 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3307 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3308 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3309 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3310 real*8 e_gcont,fprim_gcont
3311 C********************************************************************************
3318 if (dis(i,jmm).lt.dfa_cutoff) then
3319 call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3320 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3325 ci if(istrand(i,jmm).eq.1
3326 ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3329 yy1=-(dis(ipp,j)-ulhb)/dlhb
3330 y1x=rx(ipp,j)/dis(ipp,j)
3331 y1y=ry(ipp,j)/dis(ipp,j)
3332 y1z=rz(ipp,j)/dis(ipp,j)
3337 yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3338 yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3339 yy3=-pin2(i,jmm)/dshe
3340 y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3341 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3342 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3344 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3345 yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3346 yy4=-pin3(i,jmm)/dshe
3347 y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3348 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3349 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3351 yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3352 yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3353 yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3354 yy5=-pin4(i,jmm)/dshe
3355 y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3356 $ -yyy5b*rx(jm,j))*yy5
3357 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3358 $ -yyy5b*ry(jm,j))*yy5
3359 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3360 $ -yyy5b*rz(jm,j))*yy5
3373 shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3374 $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3375 shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3376 $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3377 shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3378 $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3380 ! shefx(j,11)=shefx(j,11)
3381 ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3382 ! shefy(j,11)=shefy(j,11)
3383 ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3384 ! shefz(j,11)=shefz(j,11)
3385 ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3387 yy6=-(dis(ip,j)-uldhb)/dldhb
3388 y6x=rx(ip,j)/dis(ip,j)
3389 y6y=ry(ip,j)/dis(ip,j)
3390 y6z=rz(ip,j)/dis(ip,j)
3395 yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3396 yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3397 yy8=-pina1(i,jmm)/dshe
3398 y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3399 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3400 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3402 yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3403 yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3404 yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3405 yy9=-pina2(i,jmm)/dshe
3406 y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3407 $ -yyy9b*rx(jm,j))*yy9
3408 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3409 $ -yyy9b*ry(jm,j))*yy9
3410 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3411 $ -yyy9b*rz(jm,j))*yy9
3413 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3414 yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3415 yy10=-pina4(i,jmm)/dshe
3416 y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3417 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3418 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3420 sx=y66x+y8x+y9x+y10x
3421 sy=y66y+y8y+y9y+y10y
3422 sz=y66z+y8z+y9z+y10z
3431 shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3432 $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3433 shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3434 $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3435 shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3436 $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3438 ! shefx(j,11)=shefx(j,11)
3439 ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3440 ! shefy(j,11)=shefy(j,11)
3441 ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3442 ! shefz(j,11)=shefz(j,11)
3443 ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3453 c-----------------------------------------------------------------------
3454 subroutine sheetforce12
3457 parameter(maxca=800)
3458 real*8 dfa_cutoff,dfa_cutoff_delta
3459 parameter(dfa_cutoff=15.5d0)
3460 parameter(dfa_cutoff_delta=0.5d0)
3461 cc**********************************************************************
3462 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3463 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3464 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3465 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3466 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3467 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3468 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3469 real*8 rx(maxca,maxca)
3470 real*8 ry(maxca,maxca),rz(maxca,maxca)
3471 real*8 bx(maxca),by(maxca),bz(maxca)
3472 real*8 dis(maxca,maxca)
3473 real*8 shefx(maxca,12),shefy(maxca,12)
3474 real*8 shefz(maxca,12)
3475 real*8 dp45,dm45,w_beta
3476 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3477 $ c00,s00,ulnex,dnex
3478 integer inb,nmax,iselect
3479 cc**********************************************************************
3480 common /phys1/ inb,nmax,iselect
3482 common /difvec/ rx,ry,rz
3483 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3484 $ c00,s00,ulnex,dnex
3485 common /sheconst/ dp45,dm45,w_beta
3486 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3487 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3488 common /shef/ shefx,shefy,shefz
3489 ci integer istrand(maxca,maxca)
3490 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3491 ci common /shetest/ istrand,istrand_p,istrand_m
3492 cc**********************************************************************
3494 integer j,jm,jmm,ip,i,ipp,jp
3495 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3496 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3497 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3498 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3499 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3500 real*8 e_gcont,fprim_gcont
3501 !c*************************************************************************c
3507 if (dis(i,jm).lt.dfa_cutoff) then
3508 call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3509 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3514 ci if(istrand(i,jm).eq.1
3515 ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3518 yy1=-(dis(ip,j)-ulhb)/dlhb
3519 y1x=rx(ip,j)/dis(ip,j)
3520 y1y=ry(ip,j)/dis(ip,j)
3521 y1z=rz(ip,j)/dis(ip,j)
3526 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3527 yyy3=pin1(i,jm)/(dis(ip,j)**2)
3528 yy3=-pin1(i,jm)/dshe
3529 y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3530 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3531 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3532 yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3534 yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3535 yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3536 yy4=-pin2(i,jm)/dshe
3537 y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3538 $ +yyy4b*rx(j,jp))*yy4
3539 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3540 $ +yyy4b*ry(j,jp))*yy4
3541 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3542 $ +yyy4b*rz(j,jp))*yy4
3544 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3545 yyy5=pin4(i,jm)/(dis(j,jp)**2)
3546 yy5=-pin4(i,jm)/dshe
3547 y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3548 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3549 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3562 shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3563 $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3564 shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3565 $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3566 shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3567 $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3569 ! shefx(j,12)=shefx(j,12)
3570 ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3571 ! shefy(j,12)=shefy(j,12)
3572 ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3573 ! shefz(j,12)=shefz(j,12)
3574 ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3576 yy6=-(dis(ipp,j)-uldhb)/dldhb
3577 y6x=rx(ipp,j)/dis(ipp,j)
3578 y6y=ry(ipp,j)/dis(ipp,j)
3579 y6z=rz(ipp,j)/dis(ipp,j)
3584 yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3585 yyy8=pina2(i,jm)/(dis(j,jp)**2)
3586 yy8=-pina2(i,jm)/dshe
3587 y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3588 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3589 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3591 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3592 yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3593 yy9=-pina3(i,jm)/dshe
3594 y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3595 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3596 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3598 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3599 yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3600 yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3601 yy10=-pina4(i,jm)/dshe
3602 y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3603 $ +yyy10b*rx(j,jp))*yy10
3604 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3605 $ +yyy10b*ry(j,jp))*yy10
3606 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3607 $ +yyy10b*rz(j,jp))*yy10
3609 sx=y66x+y8x+y9x+y10x
3610 sy=y66y+y8y+y9y+y10y
3611 sz=y66z+y8z+y9z+y10z
3620 shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3621 $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3622 shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3623 $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3624 shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3625 $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3636 C===============================================================================