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)
74 include 'COMMON.IOUNITS'
75 include 'COMMON.CHAIN'
77 include 'COMMON.FFIELD'
78 include 'COMMON.CONTROL'
79 include 'COMMON.SETUP'
82 C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
83 C THIS SHOULD BE MODIFIED!!
90 integer ica1, ica2,ica3,ica4,ica5
91 integer ishell, inca, itmp,iitmp
96 open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
98 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
99 print *,'Error opening dist_dfa.dat file'
102 write(iout,'(a)') 'dist_dfa.dat is opened!'
104 read(iodfa, '(a)') buffer
105 C read number of restraints
106 read(iodfa, *) IDFADIS
107 read(iodfa, *) dis_inc
109 read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
128 C READ ANGLE RESTRAINTS
130 open(iodfa, file='phi_dfa.dat',status='old',err=35)
132 35 write(iout,'(a)') 'Error opening phi_dfa.dat file'
133 print *, 'Error opening phi_dfa.dat file'
137 write(iout,'(a)') 'phi_dfa.dat is opened!'
140 read(iodfa, '(a)') buffer
141 C READ NUMBER OF RESTRAINTS
142 READ(iodfa, *) IDFAPHI
143 read(iodfa,*) phi_inc
145 read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
156 read(iodfa,*) tmp1,tmp2
170 open(iodfa, file='theta_dfa.dat',status='old',err=41)
172 41 write(iout,'(a)') 'Error opening theta_dfa.dat file'
173 print *,'Error opening theta_dfa.dat file'
176 write(iout,'(a)') 'theta_dfa.dat is opened!'
178 read(iodfa, '(a)') buffer
179 C READ NUMBER OF RESTRAINTS
180 READ(iodfa, *) IDFATHE
181 read(iodfa,*) the_inc
184 read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
195 read(iodfa,*) tmp1,tmp2
207 C END of READING ANGLE RESTRAINT!
209 C NUMBER OF NEIGHBOR CAs
210 open(iodfa,file='nei_dfa.dat',status='old',err=37)
212 37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
213 print *,'Error opening nei_dfa.dat file'
216 write(iout,'(a)') 'nei_dfa.dat is opened!'
218 read(iodfa, '(a)') buffer
219 C READ NUMBER OF RESTRAINTS
220 READ(iodfa, *) idfanei
221 read(iodfa,*) nei_inc
224 read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
233 C write(*,*) 'READ NEI:',i,j,fnei(i,j)
243 C END OF NEIGHBORING CA
246 C BETA is not parallel !
248 if (wdfa_beta.ne.0.0 .and. nfgtasks.gt.1) then
249 write (iout,*) "ERRROR dfa_beta works only for FGPROCS=1"
250 print *,"ERRROR dfa_beta works only for FGPROCS=1"
253 call int_bounds(idfadis,idfadis_start,idfadis_end)
254 call int_bounds(idfaphi,idfaphi_start,idfaphi_end)
255 call int_bounds(idfathe,idfathe_start,idfathe_end)
256 call int_bounds(idfanei,idfanei_start,idfanei_end)
257 if (me.eq.king .or. .not. out1file)
258 & write (iout,*) "DFA MPI ",
259 & "idfadis ",idfadis,idfadis_start,idfadis_end,
260 & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end,
261 & "idfathe ",idfathe,idfathe_start,idfathe_end,
262 & "idfanei ",idfanei,idfanei_start,idfanei_end
275 C READ BETA RESTRAINT
276 if (wdfa_beta.eq.0.0) return
277 open(iodfa, file='beta_dfa.dat',status='old',err=39)
279 39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
280 print *,'Error opening beta_dfa.dat file'
283 write(iout,'(a)') 'beta_dfa.dat is opened!'
285 read(iodfa,'(a)') buffer
287 read(iodfa,*) beta_inc
290 read(iodfa,*) ica1, iitmp
294 c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
299 C END OF BETA RESTRAINT
305 subroutine edfad(edfadis)
307 implicit real*8 (a-h,o-z)
309 include 'COMMON.CHAIN'
310 include 'COMMON.DERIV'
313 double precision edfadis
314 integer i, iatm1, iatm2,idiff
315 double precision ckk, sckk,dist,texp
316 double precision jix,jiy,jiz,ep,fp,scc
321 do i=idfadis_start,idfadis_end
323 iatm1=idislis(1,i)+ishiftca
324 iatm2=idislis(2,i)+ishiftca
325 idiff = abs(iatm1-iatm2)
327 JIX=c(1,iatm2)-c(1,iatm1)
328 JIY=c(2,iatm2)-c(2,iatm1)
329 JIZ=c(3,iatm2)-c(3,iatm1)
330 DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
343 if (dtmp.ge.15.0d0) then
346 c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
347 texp = exp(-dtmp)/sckk
350 ep=ep+sccdist(i,j)*texp
351 fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
353 C write(*,'(2i8,6f12.5)') i, j, dist,
354 C & fdist(i,j), ep, fp, sccdist(i,j), scc
362 c IF(ABS(EP).lt.1.0d-20)THEN
365 c IF (ABS(FP).lt.1.0d-20) THEN
369 edfadis=edfadis+ep*dis_inc*wwdist
371 gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
372 gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
373 gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
375 gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
376 gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
377 gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
384 subroutine edfat(edfator)
386 implicit real*8 (a-h,o-z)
388 include 'COMMON.CHAIN'
389 include 'COMMON.DERIV'
394 double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
395 double precision cwidth, cwidth2
396 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
404 do i=idfaphi_start,idfaphi_end
408 iatom(iii)=iphilis(iii,i)+ishiftca
411 C ANGLE VECTOR CALCULTION
412 RIX=C(1,IATOM(2))-C(1,IATOM(1))
413 RIY=C(2,IATOM(2))-C(2,IATOM(1))
414 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
416 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
417 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
418 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
420 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
421 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
422 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
424 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
425 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
426 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
428 GIX=RIY*RIPZ-RIZ*RIPY
429 GIY=RIZ*RIPX-RIX*RIPZ
430 GIZ=RIX*RIPY-RIY*RIPX
432 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
433 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
434 GIPZ=RIPX*RIPPY-RIPY*RIPPX
436 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
437 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
438 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
440 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
441 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
442 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
444 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
445 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
446 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
448 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
449 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
450 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
451 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
453 C END OF ANGLE VECTOR CALCULTION
455 TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
456 APHI(1)=TDOT/(DGI*DRIPP)
457 TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
458 APHI(2)=TDOT/(DGIP*DRIP3)
466 DDPS1=APHI(1)-FPHI1(i,j)
467 DDPS2=APHI(2)-FPHI2(i,j)
469 DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
471 if (dtmp.ge.15.0d0) then
474 c ps_tmp = dfaexp(idint(dtmp*1000)+1)
478 ephi=ephi+sccphi(i,j)*ps_tmp
480 tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
481 tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
484 C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
485 C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
488 ephi=-ephi/scc*phi_inc*wwangle
489 tfphi1=tfphi1/scc*phi_inc*wwangle
490 tfphi2=tfphi2/scc*phi_inc*wwangle
492 IF (ABS(EPHI).LT.1d-20) THEN
495 IF (ABS(TFPHI1).LT.1d-20) THEN
498 IF (ABS(TFPHI2).LT.1d-20) THEN
502 C FORCE DIRECTION CALCULATION
507 DM1=1.0d0/(DGI*DRIPP)
509 GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
510 DM2=GIRPP/(DGI**3*DRIPP)
511 DM3=GIRPP/(DGI*DRIPP**3)
513 DM4=1.0d0/(DGIP*DRIP3)
515 GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
516 DM5=GIRP3/(DGIP**3*DRIP3)
517 DM6=GIRP3/(DGIP*DRIP3**3)
519 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
520 & +( GIZ* RIPY- GIY* RIPZ)*DM2
521 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
522 & +( GIX* RIPZ- GIZ* RIPX)*DM2
523 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
524 & +( GIY* RIPX- GIX* RIPY)*DM2
528 C SECOND ATOM BY PHI1
529 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
530 & -(CIPY*GIZ-CIPZ*GIY)*DM2
531 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
532 & -(CIPZ*GIX-CIPX*GIZ)*DM2
533 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
534 & -(CIPX*GIY-CIPY*GIX)*DM2
538 C SECOND ATOM BY PHI2
540 & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
541 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
543 & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
544 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
546 & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
547 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
549 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
550 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
551 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
552 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
553 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
554 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
560 & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
561 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
563 & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
564 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
566 & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
567 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
568 C FOURTH ATOM BY PHI1
569 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
570 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
571 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
572 C FOURTH ATOM BY PHI2
574 & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
575 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
576 & + RIP3X*DM6)*TFPHI2
578 & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
579 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
580 & + RIP3Y*DM6)*TFPHI2
582 & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
583 & -( GIPX*RIPY-RIPX*GIPY)*DM5
584 & + RIP3Z*DM6)*TFPHI2
586 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
587 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
588 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
589 C END OF FORCE DIRECTION
592 gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
593 gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
594 gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
597 enephi = enephi + ephi
598 c end of single assignment statement
600 C END OF PHI RESTRAINT
602 C START OF THETA ANGLE
603 do i=idfathe_start,idfathe_end
607 iatom(iii)=ithelis(iii,i)+ishiftca
611 C ANGLE VECTOR CALCULTION
612 RIX=C(1,IATOM(2))-C(1,IATOM(1))
613 RIY=C(2,IATOM(2))-C(2,IATOM(1))
614 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
616 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
617 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
618 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
620 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
621 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
622 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
624 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
625 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
626 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
628 GIX=RIY*RIPZ-RIZ*RIPY
629 GIY=RIZ*RIPX-RIX*RIPZ
630 GIZ=RIX*RIPY-RIY*RIPX
632 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
633 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
634 GIPZ=RIPX*RIPPY-RIPY*RIPPX
636 GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
637 GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
638 GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
640 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
641 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
642 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
644 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
645 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
646 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
648 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
649 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
650 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
652 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
653 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
654 DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
655 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
656 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
657 C END OF ANGLE VECTOR CALCULTION
659 TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
660 ATHE(1)=TDOT/(DGI*DGIP)
661 TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
662 ATHE(2)=TDOT/(DGIP*DGIPP)
671 ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
672 ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
673 dtmp= (ddth1**2+ddth2**2)/cwidth2
674 if ( dtmp .ge. 15.0d0) then
677 c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
681 ethe=ethe+sccthe(i,j)*th_tmp
683 tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
684 tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
686 C write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
687 C & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
690 ethe=-ethe/scc*the_inc*wwangle
691 tfthe1=tfthe1/scc*the_inc*wwangle
692 tfthe2=tfthe2/scc*the_inc*wwangle
694 IF (ABS(ETHE).LT.TENM20) THEN
697 IF (ABS(TFTHE1).LT.TENM20) THEN
700 IF (ABS(TFTHE2).LT.TENM20) THEN
709 DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
710 DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
712 DM4=1.0d0/(DGIP*DGIPP)
713 DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
714 DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
716 C FIRST ATOM BY THETA1
717 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
718 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
719 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
720 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
721 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
722 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
723 C SECOND ATOM BY THETA1
724 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
725 & -(CIPY*GIZ-CIPZ*GIY)*DM2
726 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
727 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
728 & -(CIPZ*GIX-CIPX*GIZ)*DM2
729 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
730 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
731 & -(CIPX*GIY-CIPY*GIX)*DM2
732 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
733 C SECOND ATOM BY THETA2
735 & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
736 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
738 & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
739 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
741 & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
742 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
743 C THIRD ATOM BY THETA1
744 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
745 & -(GIY*RIZ-GIZ*RIY)*DM2
746 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
747 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
748 & -(GIZ*RIX-GIX*RIZ)*DM2
749 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
750 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
751 & -(GIX*RIY-GIY*RIX)*DM2
752 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
753 C THIRD ATOM BY THETA2
755 & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
756 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
757 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
759 & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
760 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
761 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
763 & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
764 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
765 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
766 C FOURTH ATOM BY THETA1
767 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
768 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
769 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
770 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
771 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
772 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
773 C FOURTH ATOM BY THETA2
775 & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
776 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
777 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
779 & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
780 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
781 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
783 & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
784 & -(GIPX*RIPY-GIPY*RIPX)*DM5
785 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
786 C FIFTH ATOM BY THETA2
787 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
788 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
789 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
790 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
791 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
792 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
793 C !! END OF FORCE DIRECTION!!!!
795 gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
796 gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
797 gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
800 enethe = enethe + ethe
803 edfator = enephi + enethe
808 subroutine edfan(edfanei)
809 C DFA neighboring CA restraint
810 implicit real*8 (a-h,o-z)
812 include 'COMMON.CHAIN'
813 include 'COMMON.DERIV'
817 integer kshnum, n1atom
819 double precision enenei,tmp_n
820 double precision pai,hpai
821 double precision jix,jiy,jiz,ndiff,snorm_nei
822 double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
823 double precision dr,dr2,half,ntmp,dtmp
825 parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
826 parameter(pai=3.14159265358979323846D0)
827 parameter(hpai=1.5707963267948966D0)
828 parameter(snorm_nei=0.886226925452758D0)
834 c print*, 's1:', s1(:)
835 c print*, 's2:', s2(:)
837 do i=idfanei_start,idfanei_end
840 n1atom=ineilis(i)+ishiftca
841 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
854 do j = ishiftca+1, ilastca
856 if (n1atom.eq.j) cycle
858 jix=c(1,j)-c(1,n1atom)
859 jiy=c(2,j)-c(2,n1atom)
860 jiz=c(3,j)-c(3,n1atom)
861 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
863 c write(*,*) n1atom, j, dist
866 if (dist.lt.s1(kshnum).and.
867 & dist.gt.s2(kshnum-1)) then
871 c write(*,*) 'case1:',tmp_n
880 elseif(dist.ge.s1(kshnum).and.
881 & dist.le.s2(kshnum)) then
883 dnei=(dist-s1(kshnum))/dr2*pai
884 tmp_n=tmp_n + half*(1+cos(dnei))
885 c write(*,*) 'case2:',tmp_n
886 ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
896 elseif(dist.ge.s1(kshnum-1).and.
897 & dist.le.s2(kshnum-1)) then
898 dnei=(dist-s1(kshnum-1))/dr2*pai
899 tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
900 c write(*,*) 'case3:',tmp_n
901 ftmp = hpai*sin(dnei)/dr2/dist
913 elseif(kshnum.eq.1) then
915 if(dist.lt.s1(kshnum))then
918 c write(*,*) 'case4:',tmp_n
926 elseif(dist.ge.s1(kshnum).and.
927 & dist.le.s2(kshnum))then
929 dnei=(dist-s1(kshnum))/dr2*pai
930 tmp_n=tmp_n + half*(1+cos(dnei))
931 c write(*,*) 'case5:',tmp_n
932 ftmp = -hpai*sin(dnei)/dr2/dist
953 ndiff = tmp_n-fnei(i,imin)
956 if (dtmp.ge.15.0d0) then
959 c ntmp = dfaexp( idint(dtmp*1000) + 1 )
963 enei=enei+sccnei(i,imin)*ntmp
965 & sccnei(i,imin)*ntmp*ndiff*2.0d0
966 scc=scc+sccnei(i,imin)
968 c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
969 c & fnei(i,imin),sccnei(i,imin),enei,scc
972 enei=-enei/scc*snorm_nei*nei_inc*wwnei
973 tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
975 c if (abs(enei).lt.1.0d-20)then
978 c if (abs(tmp_fnei).lt.1.0d-20) then
987 do j=ishiftca+1,ilastca
988 t2dx(j)=t2dx(j)*tmp_fnei
989 t2dy(j)=t2dy(j)*tmp_fnei
990 t2dz(j)=t2dz(j)*tmp_fnei
993 gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
994 gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
995 gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
997 do j=ishiftca+1,ilastca
998 gdfan(1,j)=gdfan(1,j)+t2dx(j)
999 gdfan(2,j)=gdfan(2,j)+t2dy(j)
1000 gdfan(3,j)=gdfan(3,j)+t2dz(j)
1002 c energy calculation
1013 subroutine edfab(edfabeta)
1015 implicit real*8 (a-h,o-z)
1017 include 'DIMENSIONS'
1018 include 'COMMON.CHAIN'
1019 include 'COMMON.DERIV'
1020 include 'COMMON.DFA'
1023 parameter(PAI=3.14159265358979323846D0)
1024 parameter (maxca=800)
1026 real*8 bx(maxres),by(maxres),bz(maxres)
1027 real*8 vbet(maxres,maxres)
1028 real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
1029 real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
1030 real*8 vbeta,vbetp,vbetm
1031 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
1032 & c00,s00,ulnex,dnex
1033 real*8 dp45,dm45,w_beta
1035 real*8 cph(maxca),cth(maxca)
1036 real*8 atx(maxca),aty(maxca),atz(maxca)
1037 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1038 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1039 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1041 real*8 astx(maxca),asty(maxca),astz(maxca)
1042 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1043 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1044 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1046 real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
1047 & astxnum(maxca),astynum(maxca),astznum(maxca),
1048 & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
1049 & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
1050 & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
1051 & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
1052 & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
1053 & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
1054 & cth_orig(maxca),sth_orig(maxca)
1056 common /sheca/ bx,by,bz
1057 common /shee/ vbeta,vbet,vbetp,vbetm
1058 common /shetf/ shetfx,shetfy,shetfz
1059 common /shef/ shefx, shefy, shefz
1060 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
1061 & c00,s00,ulnex,dnex
1062 common /sheconst/ dp45,dm45,w_beta
1064 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1065 $ atmmz,atm3x,atm3y,atm3z
1066 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1067 $ astmmz,astm3x,astm3y,astm3z
1069 common /coscos/ cph,cth
1072 C End of sheet variables
1075 double precision enebet
1078 c bx=0.0d0;by=0.0d0;bz=0.0d0
1079 c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1083 do i=ishiftca+1,ilastca
1084 bx(i-ishiftca)=c(1,i)
1085 by(i-ishiftca)=c(2,i)
1086 bz(i-ishiftca)=c(3,i)
1089 c do i=1,ilastca-ishiftca
1090 c read(99,*) bx(i),by(i),bz(i)
1098 ULNEX=COS(60.0D0/180.0D0*PAI)
1105 C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1106 S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1112 C END OF INITIALIZATION
1114 nca=ilastca-ishiftca
1116 call angvectors(nca)
1117 call sheetforce(nca,wshet)
1119 c end of sheet energy and force
1122 shetfx(j)=shetfx(j)*beta_inc
1123 shetfy(j)=shetfy(j)*beta_inc
1124 shetfz(j)=shetfz(j)*beta_inc
1125 c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1128 vbeta=vbeta*beta_inc
1133 gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1134 gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1135 gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1140 write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1154 call angvectors(nca)
1156 call angvectors(nca)
1157 atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1158 astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1160 atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1161 astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1164 atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1165 astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1168 atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1169 astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1173 call angvectors(nca)
1175 call angvectors(nca)
1177 atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1178 astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1180 atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1181 astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1184 atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1185 astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1188 atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1189 astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1193 call angvectors(nca)
1195 call angvectors(nca)
1198 atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1199 astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1201 atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1202 astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1205 atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1206 astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1209 atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1210 astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1216 write (*,'(2i5,a2,6f10.5)')
1217 & i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1218 & astxnum(i),astx(i),astxnum(i)/astx(i),
1219 & i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1220 & astynum(i),asty(i),astynum(i)/asty(i),
1221 & i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1222 & astznum(i),astz(i),astznum(i)/astz(i),
1223 & i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1224 & astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1225 & i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1226 & astmynum(i),astmy(i),astmynum(i)/astmy(i),
1227 & i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1228 & astmznum(i),astmz(i),astmznum(i)/astmz(i),
1229 & i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1230 & astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1231 & i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1232 & astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1233 & i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1234 & astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1235 & i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1236 & astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1237 & i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1238 & astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1239 & i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1240 & astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1241 & i,0," ",cth_orig(i),sth_orig(i)
1251 call angvectors(nca)
1252 call sheetforce(nca,wshet)
1253 vbeta1=vbeta*beta_inc
1255 call angvectors(nca)
1256 call sheetforce(nca,wshet)
1257 vbeta2=vbeta*beta_inc
1258 gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1262 call angvectors(nca)
1263 call sheetforce(nca,wshet)
1264 vbeta1=vbeta*beta_inc
1266 call angvectors(nca)
1267 call sheetforce(nca,wshet)
1268 vbeta2=vbeta*beta_inc
1269 gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1273 call angvectors(nca)
1274 call sheetforce(nca,wshet)
1275 vbeta1=vbeta*beta_inc
1277 call angvectors(nca)
1278 call sheetforce(nca,wshet)
1279 vbeta2=vbeta*beta_inc
1280 gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1287 call angvectors(nca)
1288 call sheetforce(nca,wshet)
1290 shetfx(j)=shetfx(j)*beta_inc
1291 shetfy(j)=shetfy(j)*beta_inc
1292 shetfz(j)=shetfz(j)*beta_inc
1296 write(*,*) 'xyz analytical and numerical gradient'
1298 write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1299 & ,(-gdfab(i,j),i=1,3)
1303 write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1304 & shetfy(j)/gdfab(2,j),
1305 & shetfz(j)/gdfab(3,j)
1313 C-------------------------------------------------------------------------------
1314 subroutine angvectors(nca)
1315 c implicit real*4(a-h,o-z)
1319 parameter(maxca=800)
1321 parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1323 real*8 bx(maxca),by(maxca),bz(maxca)
1324 real*8 dis(maxca,maxca)
1325 real*8 apx(maxca),apy(maxca),apz(maxca)
1326 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
1327 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
1328 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
1329 real*8 atx(maxca),aty(maxca),atz(maxca)
1330 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1331 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1332 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1333 real*8 astx(maxca),asty(maxca),astz(maxca)
1334 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1335 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1336 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1338 real*8 cph(maxca),cth(maxca)
1341 integer i, ip, ipp, ip3, j
1342 real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1343 real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1344 real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1345 real*8 cix, ciy, ciz, cipx, cipy, cipz
1346 real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1347 real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
1348 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1349 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1350 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1351 real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1352 real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1353 real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1354 real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1355 real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1356 real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1357 real*8 grpp,gx,gy,gz
1358 real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1359 real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1360 integer inb,nmax,iselect
1362 common /sheca/ bx,by,bz
1363 common /difvec/ rx, ry, rz
1364 common /ulang/ ulcos
1365 common /phys1/ inb,nmax,iselect
1368 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1369 & apmmz,apm3x,apm3y,apm3z
1370 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1371 & atmmz,atm3x,atm3y,atm3z
1372 common /coscos/ cph,cth
1373 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1374 & astmmz,astm3x,astm3y,astm3z
1376 C-------------------------------------------------------------------------------
1377 c write(*,*) 'inside angvectors'
1382 cph=zero; cth=zero; sth=zero
1383 apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1384 apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1385 atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1386 atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1387 astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1388 astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1391 C r[x,y,z] calc and distance calculation
1392 rx=zero;ry=zero;rz=zero
1399 dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1400 c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1401 c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1402 c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1403 c write(*,*) 'dis(i,j):',i,j,dis(i,j)
1406 c end of r[x,y,z] calc
1412 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1413 ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1414 $ +rz(i,ip)*rz(ip,ipp)
1415 ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1418 c end of virtual bond angle
1419 c write(*,*) 'inside angvectors1'
1431 rippx=bx(ip3)-bx(ipp)
1432 rippy=by(ip3)-by(ipp)
1433 rippz=bz(ip3)-bz(ipp)
1435 gx=riy*ripz-riz*ripy
1436 gy=riz*ripx-rix*ripz
1437 gz=rix*ripy-riy*ripx
1438 gpx=ripy*rippz-ripz*rippy
1439 gpy=ripz*rippx-ripx*rippz
1440 gpz=ripx*rippy-ripy*rippx
1441 gpcrp_x=gpy*ripz-gpz*ripy
1442 gpcrp_y=gpz*ripx-gpx*ripz
1443 gpcrp_z=gpx*ripy-gpy*ripx
1444 d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1445 gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1446 & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1452 gmx=rimy*riz-rimz*riy
1453 gmy=rimz*rix-rimx*riz
1454 gmz=rimx*riy-rimy*rix
1455 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1457 ggm=gmx*gx+gmy*gy+gmz*gz
1458 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1464 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1466 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1467 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1469 c write(*,*) 'inside angvectors2'
1471 rimmx=bx(i-1)-bx(i-2)
1472 rimmy=by(i-1)-by(i-2)
1473 rimmz=bz(i-1)-bz(i-2)
1475 gmmx=rimmy*rimz-rimmz*rimy
1476 gmmy=rimmz*rimx-rimmx*rimz
1477 gmmz=rimmx*rimy-rimmy*rimx
1478 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1480 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1481 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1482 gmcrim_x=gmy*rimz-gmz*rimy
1483 gmcrim_y=gmz*rimx-gmx*rimz
1484 gmcrim_z=gmx*rimy-gmy*rimx
1485 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1486 d_gmcrim3=d_gmcrim**3
1487 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1488 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1492 rim3x=bx(i-2)-bx(i-3)
1493 rim3y=by(i-2)-by(i-3)
1494 rim3z=bz(i-2)-bz(i-3)
1495 g3x=rim3y*rimmz-rim3z*rimmy
1496 g3y=rim3z*rimmx-rim3x*rimmz
1497 g3z=rim3x*rimmy-rim3y*rimmx
1498 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1499 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1500 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1501 cc**********************************************************************
1502 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1503 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1504 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1505 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1506 d_gmmcrimm3=d_gmmcrimm**3
1507 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1508 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1515 dg=sqrt(gx**2+gy**2+gz**2)
1516 dgp=sqrt(gpx**2+gpy**2+gpz**2)
1519 ggp=gx*gpx+gy*gpy+gz*gpz
1520 grpp=gx*rippx+gy*rippy+gz*rippz
1522 if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1523 & .and.d_gpcrp.gt.0.0D0) then
1524 cph(i)=grpp/dg/dripp
1526 sth(i)=gpcrp__g/d_gpcrp/dg
1534 c write(*,*) 'inside angvectors3'
1536 if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1537 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1540 d12=1.0D0/(dg*dripp)
1541 d13=grpp/(dg3*dripp)
1542 sd10=1.0D0/(d_gpcrp*dg)
1543 sd11=gpcrp__g/(d_gpcrp*dg3)
1553 atx(i)=(ripz*gpy-ripy*gpz)*d10
1554 & -(gy*ripz-gz*ripy)*d11
1555 aty(i)=(ripx*gpz-ripz*gpx)*d10
1556 & -(gz*ripx-gx*ripz)*d11
1557 atz(i)=(ripy*gpx-ripx*gpy)*d10
1558 & -(gx*ripy-gy*ripx)*d11
1559 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1560 & +ripy*gpy*ripx-gpx*ripz**2)
1561 & -sd11*(gy*ripz-gz*ripy)
1562 asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1563 & -gpy*ripx**2+gpz*ripy*ripz)
1564 & -sd11*(-gx*ripz+gz*ripx)
1565 astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1566 & -gpz*ripy**2+ripz*gpx*ripx)
1567 & -sd11*(gx*ripy-gy*ripx)
1568 apx(i)=(ripz*rippy-ripy*rippz)*d12
1569 & -(gy*ripz-gz*ripy)*d13
1570 apy(i)=(ripx*rippz-ripz*rippx)*d12
1571 & -(gz*ripx-gx*ripz)*d13
1572 apz(i)=(ripy*rippx-ripx*rippy)*d12
1573 & -(gx*ripy-gy*ripx)*d13
1585 if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1586 & .and.d_gcr3.gt.0.0D0) then
1590 d23=1.0D0/(dgm*drip)
1591 d24=gmrp/(dgm3*drip)
1592 sd20=1.0D0/(d_gcr*dgm)
1593 sd21=gcr__gm/(d_gcr3*dgm)
1594 sd22=gcr__gm/(d_gcr*dgm3)
1605 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1606 & -(ciy*gmz-ciz*gmy)*d21
1607 & +(ripy*gz-ripz*gy)*d22
1608 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1609 & -(ciz*gmx-cix*gmz)*d21
1610 & +(ripz*gx-ripx*gz)*d22
1611 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1612 & -(cix*gmy-ciy*gmx)*d21
1613 & +(ripx*gy-ripy*gx)*d22
1614 cc**********************************************************************
1615 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1616 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1617 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1618 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1619 & +gcr_z*(-ripz*rix+gy))
1620 & -sd22*(-gmy*ciz+gmz*ciy)
1622 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1623 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1625 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1626 & -gcr_z*(ripz*riy+gx))
1627 & -sd22*(gmx*ciz-gmz*cix)
1629 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1630 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1632 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1633 & +gcr_z*(ripy*riy+ripx*rix))
1634 & -sd22*(-gmx*ciy+gmy*cix)
1635 cc**********************************************************************
1636 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1637 & -(ciy*gmz-ciz*gmy)*d24
1638 apmy(i)=(ciz*ripx-ripz*cix)*d23
1639 & -(ciz*gmx-cix*gmz)*d24
1640 apmz(i)=(cix*ripy-ripx*ciy)*d23
1641 & -(cix*gmy-ciy*gmx)*d24
1645 if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1646 & .and.d_gmcrim3.gt.0.0D0) then
1647 d30=1.0D0/(dgm*dgmm)
1648 d31=gmmgm/(dgm3*dgmm)
1649 d32=gmmgm/(dgm*dgmm3)
1650 d33=1.0D0/(dgmm*dri)
1651 d34=gmmr/(dgmm3*dri)
1652 d35=gmmr/(dgmm*dri3)
1653 sd30=1.0D0/(d_gmcrim*dgmm)
1654 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1655 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1668 c write(*,*) 'inside angvectors4'
1670 cc**********************************************************************
1671 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1672 & -(ciy*gmz-ciz*gmy)*d31
1673 & -(gmmy*rimmz-gmmz*rimmy)*d32
1674 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1675 & -(ciz*gmx-cix*gmz)*d31
1676 & -(gmmz*rimmx-gmmx*rimmz)*d32
1677 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1678 & -(cix*gmy-ciy*gmx)*d31
1679 & -(gmmx*rimmy-gmmy*rimmx)*d32
1680 cc**********************************************************************
1681 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1682 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1683 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1684 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1685 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1686 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1687 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1689 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1690 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1691 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1692 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1693 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1694 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1695 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1697 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1698 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1699 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1700 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1701 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1702 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1703 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1704 c**********************************************************************
1705 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1706 & -(gmmy*rimmz-gmmz*rimmy)*d34
1708 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1709 & -(gmmz*rimmx-gmmx*rimmz)*d34
1711 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1712 & -(gmmx*rimmy-gmmy*rimmx)*d34
1717 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1718 & .and.drim3.gt.0.0D0
1719 & .and.d_gmmcrimm3.gt.0.0D0) then
1720 d40=1.0D0/(dg30*dgmm)
1721 d41=g3gmm/(dg30*dgmm3)
1722 d42=1.0D0/(dg30*drim)
1723 d43=g3rim_/(dg30*drim3)
1724 sd40=1.0D0/(dg30*d_gmmcrimm)
1725 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1734 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1735 & -(gmmy*rimmz-gmmz*rimmy)*d41
1736 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1737 & -(gmmz*rimmx-gmmx*rimmz)*d41
1738 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1739 & -(gmmx*rimmy-gmmy*rimmx)*d41
1740 cc**********************************************************************
1741 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1742 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1743 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1744 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1746 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1747 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1748 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1749 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1751 c & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1753 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1754 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1755 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1756 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1757 c**********************************************************************
1758 apm3x(i)=g3x*d42-rimx*d43
1759 apm3y(i)=g3y*d42-rimy*d43
1760 apm3z(i)=g3z*d42-rimz*d43
1763 c*******************************************************************************
1765 c write(*,*) 'inside angvectors5'
1772 rimmx=bx(i-1)-bx(i-2)
1773 rimmy=by(i-1)-by(i-2)
1774 rimmz=bz(i-1)-bz(i-2)
1775 rim3x=bx(i-2)-bx(i-3)
1776 rim3y=by(i-2)-by(i-3)
1777 rim3z=bz(i-2)-bz(i-3)
1778 gmmx=rimmy*rimz-rimmz*rimy
1779 gmmy=rimmz*rimx-rimmx*rimz
1780 gmmz=rimmx*rimy-rimmy*rimx
1781 g3x=rim3y*rimmz-rim3z*rimmy
1782 g3y=rim3z*rimmx-rim3x*rimmz
1783 g3z=rim3x*rimmy-rim3y*rimmx
1785 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1786 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1787 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1792 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1793 cc**********************************************************************
1794 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1795 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1796 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1797 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1798 d_gmmcrimm3=d_gmmcrimm**3
1799 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1800 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1802 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1803 & .and.drim3.gt.0.0D0
1804 & .and.d_gmmcrimm3.gt.0.0D0) then
1805 d40=1.0D0/(dg30*dgmm)
1806 d41=g3gmm/(dg30*dgmm3)
1807 d42=1.0D0/(dg30*drim)
1808 d43=g3rim_/(dg30*drim3)
1809 sd40=1.0D0/(dg30*d_gmmcrimm)
1810 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1819 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1820 & -(gmmy*rimmz-gmmz*rimmy)*d41
1821 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1822 & -(gmmz*rimmx-gmmx*rimmz)*d41
1823 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1824 & -(gmmx*rimmy-gmmy*rimmx)*d41
1825 cc**********************************************************************
1826 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1827 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1828 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1829 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1831 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1832 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1833 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1834 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1836 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1837 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1838 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1839 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1840 cc**********************************************************************
1841 apm3x(i)=g3x*d42-rimx*d43
1842 apm3y(i)=g3y*d42-rimy*d43
1843 apm3z(i)=g3z*d42-rimz*d43
1853 gmx=rimy*riz-rimz*riy
1854 gmy=rimz*rix-rimx*riz
1855 gmz=rimx*riy-rimy*rix
1856 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1860 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1861 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1862 gmcrim_x=gmy*rimz-gmz*rimy
1863 gmcrim_y=gmz*rimx-gmx*rimz
1864 gmcrim_z=gmx*rimy-gmy*rimx
1865 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1866 d_gmcrim3=d_gmcrim**3
1867 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1868 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1870 if(dgm3.gt.0.0D0.and.
1871 & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1872 & .and.d_gmcrim3.gt.0.0D0) then
1873 d30=1.0D0/(dgm*dgmm)
1874 d31=gmmgm/(dgm3*dgmm)
1875 d32=gmmgm/(dgm*dgmm3)
1876 d33=1.0D0/(dgmm*dri)
1877 d34=gmmr/(dgmm3*dri)
1878 d35=gmmr/(dgmm*dri3)
1879 sd30=1.0D0/(d_gmcrim*dgmm)
1880 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1881 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1894 cc**********************************************************************
1895 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1896 & -(ciy*gmz-ciz*gmy)*d31
1897 & -(gmmy*rimmz-gmmz*rimmy)*d32
1898 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1899 & -(ciz*gmx-cix*gmz)*d31
1900 & -(gmmz*rimmx-gmmx*rimmz)*d32
1901 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1902 & -(cix*gmy-ciy*gmx)*d31
1903 & -(gmmx*rimmy-gmmy*rimmx)*d32
1904 cc**********************************************************************
1905 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1906 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1907 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1908 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1909 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1910 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1911 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1913 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1914 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1915 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1916 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1917 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1918 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1919 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1921 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1922 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1923 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1924 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1925 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1926 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1927 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1928 cc**********************************************************************
1929 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1930 & -(gmmy*rimmz-gmmz*rimmy)*d34
1932 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1933 & -(gmmz*rimmx-gmmx*rimmz)*d34
1935 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1936 & -(gmmx*rimmy-gmmy*rimmx)*d34
1940 c write(*,*) 'inside angvectors6'
1950 gx=riy*ripz-riz*ripy
1951 gy=riz*ripx-rix*ripz
1952 gz=rix*ripy-riy*ripx
1953 ggm=gmx*gx+gmy*gy+gmz*gz
1954 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1955 dg=sqrt(gx**2+gy**2+gz**2)
1961 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1963 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1964 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1965 if(dgm3.gt.0.0D0.and.
1966 & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1971 d23=1.0D0/(dgm*drip)
1972 d24=gmrp/(dgm3*drip)
1973 sd20=1.0D0/(d_gcr*dgm)
1974 sd21=gcr__gm/(d_gcr3*dgm)
1975 sd22=gcr__gm/(d_gcr*dgm3)
1986 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1987 & -(ciy*gmz-ciz*gmy)*d21
1988 & +(ripy*gz-ripz*gy)*d22
1989 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1990 & -(ciz*gmx-cix*gmz)*d21
1991 & +(ripz*gx-ripx*gz)*d22
1992 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1993 & -(cix*gmy-ciy*gmx)*d21
1994 & +(ripx*gy-ripy*gx)*d22
1995 cc**********************************************************************
1996 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1997 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1998 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1999 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
2000 & +gcr_z*(-ripz*rix+gy))
2001 & -sd22*(-gmy*ciz+gmz*ciy)
2003 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
2004 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
2006 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
2007 & -gcr_z*(ripz*riy+gx))
2008 & -sd22*(gmx*ciz-gmz*cix)
2010 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
2011 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
2013 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
2014 & +gcr_z*(ripy*riy+ripx*rix))
2015 & -sd22*(-gmx*ciy+gmy*cix)
2016 cc**********************************************************************
2018 apmx(i)=(ciy*ripz-ripy*ciz)*d23
2019 & -(ciy*gmz-ciz*gmy)*d24
2020 apmy(i)=(ciz*ripx-ripz*cix)*d23
2021 & -(ciz*gmx-cix*gmz)*d24
2022 apmz(i)=(cix*ripy-ripx*ciy)*d23
2023 & -(cix*gmy-ciy*gmx)*d24
2031 c-------------------------------------------------------------------------------
2032 C---------------------------------------------------------------------------------
2033 subroutine sheetforce(nca,wshet)
2036 c this should be matched with dfa.fcm
2038 parameter(maxca=800)
2039 cc**********************************************************************
2042 integer inb,nmax,iselect
2044 c real*8 dfaexp(15001)
2046 real*8 vbeta,vbetp,vbetm
2047 real*8 shefx(maxca,12)
2048 real*8 shefy(maxca,12),shefz(maxca,12)
2049 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2050 real*8 vbet(maxca,maxca)
2051 real*8 wshet(maxca,maxca)
2052 real*8 bx(maxca),by(maxca),bz(maxca)
2054 common /sheca/ bx,by,bz
2055 common /phys1/ inb,nmax,iselect
2056 common /shef/ shefx,shefy,shefz
2057 common /shee/ vbeta,vbet,vbetp,vbetm
2058 common /shetf/ shetfx,shetfy,shetfz
2075 call sheetene(nca,wshet)
2078 887 format(a,1x,i6,3x,f12.8)
2079 888 format(a,1x,i4,1x,i4,3x,f12.8)
2080 889 format(a,1x,i4,3x,f12.8)
2081 !write(2,*) 'coord : '
2083 !write(2,887) 'bx:',i,bx(i)
2084 !write(2,887) 'by:',i,by(i)
2085 !write(2,887) 'bz:',i,bz(i)
2087 !write(2,*) 'After sheetforce1'
2090 !write(2,888) 'shefx :',i,k,shefx(i,k)
2091 !write(2,888) 'shefy :',i,k,shefy(i,k)
2092 !write(2,888) 'shefz :',i,k,shefz(i,k)
2098 !write(2,*) 'After sheetforce5'
2101 !write(2,888) 'shefx :',i,k,shefx(i,k)
2102 !write(2,888) 'shefy :',i,k,shefy(i,k)
2103 !write(2,888) 'shefz :',i,k,shefz(i,k)
2109 !write(2,*) 'After sheetforce6'
2112 !write(2,888) 'shefx :',i,k,shefx(i,k)
2113 !write(2,888) 'shefy :',i,k,shefy(i,k)
2114 !write(2,888) 'shefz :',i,k,shefz(i,k)
2120 !write(2,*) 'After sheetforce11'
2123 !write(2,888) 'shefx :',i,k,shefx(i,k)
2124 !write(2,888) 'shefy :',i,k,shefy(i,k)
2125 !write(2,888) 'shefz :',i,k,shefz(i,k)
2131 !write(2,*) 'After sheetforce12'
2134 !write(2,888) 'shefx :',i,k,shefx(i,k)
2135 !write(2,888) 'shefy :',i,k,shefy(i,k)
2136 !write(2,888) 'shefz :',i,k,shefz(i,k)
2142 shetfx(i)=shetfx(i)+shefx(i,k)
2143 shetfy(i)=shetfy(i)+shefy(i,k)
2144 shetfz(i)=shetfz(i)+shefz(i,k)
2147 !write(2,*) 'Beta Finished'
2149 !write(2,889) 'shetfx : ',i,shetfx(i)
2150 !write(2,889) 'shetfy : ',i,shetfy(i)
2151 !write(2,889) 'shetfz : ',i,shetfz(i)
2157 c-------------------------------------------------------------------------------
2158 subroutine sheetene(nca,wshet)
2161 parameter(maxca=800)
2162 real*8 dfa_cutoff,dfa_cutoff_delta
2163 parameter(dfa_cutoff=15.5d0)
2164 parameter(dfa_cutoff_delta=0.5d0)
2165 cc******************************************************************************
2167 c real*8 dfaexp(15001)
2168 real*8 dtmp1, dtmp2, dtmp3
2170 real*8 vbet(maxca,maxca)
2171 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2172 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2173 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2174 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2175 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2176 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2177 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2178 real*8 cph(maxca),cth(maxca)
2179 real*8 rx(maxca,maxca)
2180 real*8 ry(maxca,maxca),rz(maxca,maxca)
2181 real*8 bx(maxca),by(maxca),bz(maxca)
2182 real*8 dis(maxca,maxca)
2184 cc**********************************************************************
2185 real*8 astx(maxca),asty(maxca),astz(maxca)
2186 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2187 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2188 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2190 real*8 wshet(maxca,maxca)
2191 real*8 dp45, dm45, w_beta
2192 real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2194 integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2196 real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2198 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2199 common /shetf/ shetfx,shetfy,shetfz
2201 common /sheca/ bx,by,bz
2202 common /phys1/ inb,nmax,iselect
2204 common /difvec/ rx,ry,rz
2205 common /coscos/ cph,cth
2206 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2207 & c00,s00,ulnex,dnex
2208 common /sheconst/ dp45,dm45,w_beta
2209 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2210 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2211 common /shee/ vbeta,vbet,vbetp,vbetm
2212 common /ulang/ ulcos
2213 cc**********************************************************************
2214 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2215 & astmmz,astm3x,astm3y,astm3z
2218 real*8 r_pair_mat(maxca,maxca)
2219 real*8 e_gcont,fprim_gcont,de_gcont
2220 ci integer istrand(maxca,maxca)
2221 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2222 ci common /shetest/ istrand,istrand_p,istrand_m
2223 common /beta_p/ r_pair_mat
2224 C-------------------------------------------------------------------------------
2228 r_pair_mat(i,j)=wshet(i,j)
2229 c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2241 if (dis(i,j).lt.dfa_cutoff) then
2242 call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2243 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2250 cc**********************************************************************
2251 y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2252 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2254 y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2255 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2268 ci if(istrand(i,j).eq.1) then
2270 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2271 yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2274 pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2275 $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2276 pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2277 $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2278 pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2279 $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2280 pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2281 $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2283 yshe1=pin1(i,j)**2+pin2(i,j)**2
2284 yshe1=-0.5d0*yshe1/dshe
2285 yshe2=pin3(i,j)**2+pin4(i,j)**2
2286 yshe2=-0.5d0*yshe2/dshe
2288 ci if((yshe1+yshe2).ge.-4) then
2295 C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2296 C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2297 C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2298 C write(*,*) 'dis(i,j):',i,j,dis(i,j)
2299 C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2300 C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2301 C write(*,*) 'pin1:',pin1(i,j)
2302 C write(*,*) 'pin2:',pin2(i,j)
2303 C write(*,*) 'pin3:',pin3(i,j)
2304 C write(*,*) 'pin4:',pin4(i,j)
2307 C write(*,*) 'yy1:',yy1
2308 C write(*,*) 'yy2:',yy2
2309 C write(*,*) 'yshe1:',yshe1
2310 C write(*,*) 'yshe2:',yshe2
2313 ci if (istrand_p(i,j).eq.1) then
2321 dtmp3 = y+yy1+yy2+yshe1+yshe2
2323 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2324 C write(*,*)'2', y,yy1,yy2
2325 C write(*,*)'3', yshe1,yshe2
2327 cc if (dtmp3.le.-35.0d0) then
2328 c vbetap(i,j)=-dp45*exp(dtmp3)
2329 cc vbetap(i,j)=0.0d0
2331 c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2332 vbetap(i,j)=-dp45*exp(dtmp3)
2335 cc if (dtmp1.le.-35.0d0) then
2336 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2337 cc vbetap1(i,j)=0.0d0
2339 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2340 c $ *dfaexp(idint(-dtmp1*1000)+1)
2341 vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2344 cc if (dtmp2.le.-35.0d0) then
2345 C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2346 cc vbetap2(i,j)=0.0d0
2348 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2349 c $ *dfaexp(idint(-dtmp2*1000)+1)
2350 vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2353 c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2354 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2355 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2357 ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2358 ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2360 ci elseif (istrand_p(i,j).eq.0)then
2366 yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2367 yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2369 pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2370 $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2371 pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2372 $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2373 pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2374 $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2375 pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2376 $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2378 yshe1=pina1(i,j)**2+pina2(i,j)**2
2379 yshe1=-0.5d0*yshe1/dshe
2380 yshe2=pina3(i,j)**2+pina4(i,j)**2
2381 yshe2=-0.5d0*yshe2/dshe
2383 ci if((yshe1+yshe2).ge.-4) then
2390 C write(*,*) 'pina1:',pina1(i,j)
2391 C write(*,*) 'pina2:',pina2(i,j)
2392 C write(*,*) 'pina3:',pina3(i,j)
2393 C write(*,*) 'pina4:',pina4(i,j)
2394 C write(*,*) 'yshe1:',yshe1
2395 C write(*,*) 'yshe2:',yshe2
2396 C write(*,*) 'dshe:',dshe
2398 ci if (istrand_m(i,j).eq.1) then
2405 dtmp3=y+yy1+yy2+yshe1+yshe2
2409 cc if(dtmp3 .le. -35.0d0) then
2410 c vbetam(i,j)=-dm45*exp(dtmp3)
2411 cc vbetam(i,j)=0.0d0
2413 c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2414 vbetam(i,j)=-dm45*exp(dtmp3)
2417 cc if(dtmp1 .le. -35.0d0) then
2418 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2419 cc vbetam1(i,j)=0.0d0
2421 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2422 c $ *dfaexp(idint(-dtmp1*1000)+1)
2423 vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2426 cc if(dtmp2.le.-35.0d0) then
2427 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2428 cc vbetam2(i,j)=0.0d0
2430 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2431 c $ *dfaexp(idint(-dtmp2*1000)+1)
2432 vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2435 ci elseif (istrand_m(i,j).eq.0)then
2442 c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2443 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2444 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2446 ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2447 ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2449 uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2450 uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2452 c write(*,*) 'uup,uum:', uup, uum
2454 c uup=vbetap1(i,j)+vbetap2(i,j)
2455 c uum=vbetam1(i,j)+vbetam2(i,j)
2460 vbeta=vbeta+vbet(i,j)*e_gcont
2463 if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
2464 c gradient correction from gcont
2465 de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
2466 shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
2467 shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
2468 shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
2470 shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
2471 shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
2472 shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
2474 c energy correction from gcont
2475 vbet(i,j)=vbet(i,j)*e_gcont
2476 vbetap(i,j)=vbetap(i,j)*e_gcont
2477 vbetap1(i,j)=vbetap1(i,j)*e_gcont
2478 vbetap2(i,j)=vbetap2(i,j)*e_gcont
2479 vbetam(i,j)=vbetam(i,j)*e_gcont
2480 vbetam1(i,j)=vbetam1(i,j)*e_gcont
2481 vbetam2(i,j)=vbetam2(i,j)*e_gcont
2485 ci elseif(istrand(i,j).eq.0)then
2489 c write(*,*) 'uup,uum:',uup,uum
2490 c write(*,*) 'vbetap(i,j):',vbetap(i,j)
2491 c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2492 c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2493 c write(*,*) 'vbetam(i,j):',vbetam(i,j)
2494 c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2495 c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2496 c write(*,*) 'uup:',uup
2497 c write(*,*) 'uum:',uum
2498 c write(*,*) 'vbetp:',vbetp
2499 c write(*,*) 'vbetm:',vbetm
2500 c write(*,*) 'vbet(i,j):',vbet(i,j)
2517 ! write(*,*) 'I,J:', i,j
2518 ! write(*,*) 'vbetap(i,j):',vbetap(i,j)
2519 ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2520 ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2521 ! write(*,*) 'vbetam(i,j):',vbetam(i,j)
2522 ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2523 ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2524 ! write(*,*) 'vbet(i,j):',vbet(i,j)
2530 c-------------------------------------------------------------------------------
2531 subroutine sheetforce1
2534 parameter(maxca=800)
2535 real*8 dfa_cutoff,dfa_cutoff_delta
2536 parameter(dfa_cutoff=15.5d0)
2537 parameter(dfa_cutoff_delta=0.5d0)
2538 cc**********************************************************************
2539 real*8 vbet(maxca,maxca)
2540 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2541 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2542 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2543 real*8 cph(maxca),cth(maxca)
2544 real*8 rx(maxca,maxca)
2545 real*8 ry(maxca,maxca),rz(maxca,maxca)
2546 real*8 bx(maxca),by(maxca),bz(maxca)
2547 real*8 dis(maxca,maxca)
2548 real*8 shefx(maxca,12)
2549 real*8 shefy(maxca,12),shefz(maxca,12)
2550 real*8 atx(maxca),aty(maxca),atz(maxca)
2551 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2552 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2553 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2554 real*8 apx(maxca),apy(maxca),apz(maxca)
2555 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2556 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2557 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2559 real*8 astx(maxca),asty(maxca),astz(maxca)
2560 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2561 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2562 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2564 real*8 w_beta,dp45, dm45
2565 real*8 vbeta, vbetp, vbetm
2566 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2567 $ c00,s00,ulnex,dnex
2568 integer inb,nmax,iselect
2570 common /phys1/ inb,nmax,iselect
2572 common /difvec/ rx,ry,rz
2573 common /coscos/ cph,cth
2574 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2575 $ c00,s00,ulnex,dnex
2576 common /sheconst/ dp45,dm45,w_beta
2577 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2578 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2579 $ atmmz,atm3x,atm3y,atm3z
2580 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2581 $ apmmz,apm3x,apm3y,apm3z
2582 common /shef/ shefx,shefy,shefz
2583 common /shee/ vbeta,vbet,vbetp,vbetm
2584 common /ulang/ ulcos
2585 c c**********************************************************************
2586 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2587 $ astmmz,astm3x,astm3y,astm3z
2589 C--------------------------------------------------------------------------------
2591 integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2592 real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2593 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2594 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2595 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2596 real*8 e_gcont,fprim_gcont
2597 C--------------------------------------------------------------------------------
2602 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2607 cc1=(ulcos(imm)-ulnex)/dnex
2608 dmm=cc1/(dis(imm,im)*dis(im,i))
2609 dmm__=cc1*ulcos(imm)/dis(im,i)**2
2610 fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2611 fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2612 fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2616 fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2617 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2618 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2628 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2633 cc1=(ulcos(imm)-ulnex)/dnex
2634 cc2=(ulcos(im)-ulnex)/dnex
2635 dmm1=cc1/(dis(imm,im)*dis(im,i))
2636 dmm2=cc2/(dis(im,i)*dis(i,ip))
2637 dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2638 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2639 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2640 cc**********************************************************************
2641 fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2642 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2643 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2644 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2645 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2646 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2650 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2651 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2652 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2661 c3=(cth(im)*c00+sth(im)*s00-1)/dca
2666 cc2=(ulcos(im)-ulnex)/dnex
2667 cc3=(ulcos(i)-ulnex)/dnex
2668 dmm2=cc2/(dis(im,i)*dis(i,ip))
2669 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2670 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2671 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2672 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2673 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2674 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2675 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2676 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2677 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2678 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2682 fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2683 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2684 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2692 c4=(cth(i)*c00+sth(i)*s00-1)/dca
2697 cc3=(ulcos(i)-ulnex)/dnex
2698 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2699 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2700 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2701 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2702 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2706 fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2707 fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2708 fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2717 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2722 cc7=(ulcos(jmm)-ulnex)/dnex
2723 dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2724 dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2725 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2726 fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2727 fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2731 fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2732 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2733 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2742 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2747 cc7=(ulcos(jmm)-ulnex)/dnex
2748 cc8=(ulcos(jm)-ulnex)/dnex
2749 dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2750 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2751 dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2752 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2753 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2754 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2755 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2756 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2757 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2758 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2759 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2763 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2764 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2765 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2775 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2780 cc8=(ulcos(jm)-ulnex)/dnex
2781 cc9=(ulcos(j)-ulnex)/dnex
2782 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2783 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2784 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2785 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2786 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2787 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2788 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2789 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2790 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2791 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2792 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2796 fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2797 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2798 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2807 c10=(cth(j)*c00+sth(j)*s00-1)/dca
2812 cc9=(ulcos(j)-ulnex)/dnex
2813 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2814 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2815 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2816 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2817 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2821 fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2822 fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2823 fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2831 c----------------------------------------------------------------------------
2832 subroutine sheetforce5
2835 parameter(maxca=800)
2836 real*8 dfa_cutoff,dfa_cutoff_delta
2837 parameter(dfa_cutoff=15.5d0)
2838 parameter(dfa_cutoff_delta=0.5d0)
2839 cc**********************************************************************
2840 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2841 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2842 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2843 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2844 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2845 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2846 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2847 real*8 rx(maxca,maxca)
2848 real*8 ry(maxca,maxca),rz(maxca,maxca)
2849 real*8 bx(maxca),by(maxca),bz(maxca)
2850 real*8 dis(maxca,maxca)
2851 real*8 shefx(maxca,12),shefy(maxca,12)
2852 real*8 shefz(maxca,12)
2853 real*8 dp45,dm45,w_beta
2854 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2855 $ c00,s00,ulnex,dnex
2856 integer inb,nmax,iselect
2857 cc**********************************************************************
2858 common /phys1/ inb,nmax,iselect
2860 common /difvec/ rx,ry,rz
2861 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2862 $ c00,s00,ulnex,dnex
2863 common /sheconst/ dp45,dm45,w_beta
2864 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2865 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2866 common /shef/ shefx,shefy,shefz
2867 ci integer istrand(maxca,maxca)
2868 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2869 ci common /shetest/ istrand,istrand_p,istrand_m
2870 c********************************************************************************
2872 integer i,imm,im,jp,jpp,j
2873 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2874 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2875 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2876 real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2877 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2878 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2879 real*8 e_gcont,fprim_gcont
2880 c********************************************************************************
2886 if (dis(imm,j).lt.dfa_cutoff) then
2887 call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2888 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2893 ci if(istrand(imm,j).eq.1
2894 ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2897 yy1=-(dis(i,jpp)-ulhb)/dlhb
2898 y1x=rx(jpp,i)/dis(i,jpp)
2899 y1y=ry(jpp,i)/dis(i,jpp)
2900 y1z=rz(jpp,i)/dis(i,jpp)
2905 yy33=1.0D0/(dis(im,jp)*dis(im,i))
2906 yyy3=pin1(imm,j)/(dis(im,i)**2)
2907 yy3=-pin1(imm,j)/dshe
2908 y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2909 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2910 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2912 yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2913 yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2914 yyy4b=pin3(imm,j)/(dis(im,i)**2)
2915 yy4=-pin3(imm,j)/dshe
2916 y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2917 $ -yyy4b*rx(im,i))*yy4
2918 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2919 $ -yyy4b*ry(im,i))*yy4
2920 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2921 $ -yyy4b*rz(im,i))*yy4
2924 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2925 yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2926 yy5=-pin4(imm,j)/dshe
2927 y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2928 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2929 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2942 shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2943 $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2944 shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2945 $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2946 shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2947 $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2949 ! shefx(i,5)=shefx(i,5)
2950 ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2951 ! shefy(i,5)=shefy(i,5)
2952 ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2953 ! shefz(i,5)=shefz(i,5)
2954 ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2956 yy6=-(dis(i,jp)-uldhb)/dldhb
2957 y6x=rx(jp,i)/dis(i,jp)
2958 y6y=ry(jp,i)/dis(i,jp)
2959 y6z=rz(jp,i)/dis(i,jp)
2964 yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2965 yyy8=pina1(imm,j)/(dis(im,i)**2)
2966 yy8=-pina1(imm,j)/dshe
2967 y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2968 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2969 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2971 yy99=1.0D0/(dis(jp,i)*dis(im,i))
2972 yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2973 yyy9b=pina3(imm,j)/(dis(im,i)**2)
2974 yy9=-pina3(imm,j)/dshe
2975 y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2976 $ -yyy9b*rx(im,i))*yy9
2977 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2978 $ -yyy9b*ry(im,i))*yy9
2979 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2980 $ -yyy9b*rz(im,i))*yy9
2982 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2983 yyy10=pina4(imm,j)/(dis(jp,i)**2)
2984 yy10=-pina4(imm,j)/dshe
2985 y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2986 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2987 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2989 sx=y66x+y8x+y9x+y10x
2990 sy=y66y+y8y+y9y+y10y
2991 sz=y66z+y8z+y9z+y10z
3000 shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
3001 $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3002 shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
3003 $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3004 shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
3005 $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3007 ! shefx(i,5)=shefx(i,5)
3008 ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3009 ! shefy(i,5)=shefy(i,5)
3010 ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3011 ! shefz(i,5)=shefz(i,5)
3012 ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3022 c--------------------------------------------------------------------------c
3023 subroutine sheetforce6
3026 parameter(maxca=800)
3027 real*8 dfa_cutoff,dfa_cutoff_delta
3028 parameter(dfa_cutoff=15.5d0)
3029 parameter(dfa_cutoff_delta=0.5d0)
3030 cc**********************************************************************
3031 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3032 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3033 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3034 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3035 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3036 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3037 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3038 real*8 rx(maxca,maxca)
3039 real*8 ry(maxca,maxca),rz(maxca,maxca)
3040 real*8 bx(maxca),by(maxca),bz(maxca)
3041 real*8 dis(maxca,maxca)
3042 real*8 shefx(maxca,12),shefy(maxca,12)
3043 real*8 shefz(maxca,12)
3044 real*8 dp45,dm45,w_beta
3045 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3046 $ c00,s00,ulnex,dnex
3047 integer inb,nmax,iselect
3048 cc**********************************************************************
3049 common /phys1/ inb,nmax,iselect
3051 common /difvec/ rx,ry,rz
3052 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3053 $ c00,s00,ulnex,dnex
3054 common /sheconst/ dp45,dm45,w_beta
3055 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3056 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3057 common /shef/ shefx,shefy,shefz
3058 ci integer istrand(maxca,maxca)
3059 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3060 ci common /shetest/ istrand,istrand_p,istrand_m
3061 cc**********************************************************************
3063 integer i,imm,im,jp,jpp,j,ip
3064 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3065 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3066 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3067 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3068 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
3069 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
3070 real*8 e_gcont,fprim_gcont
3071 C********************************************************************************
3077 if (dis(im,j).lt.dfa_cutoff) then
3078 call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3079 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3084 ci if(istrand(im,j).eq.1
3085 ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
3088 yy1=-(dis(i,jp)-ulhb)/dlhb
3089 y1x=rx(jp,i)/dis(i,jp)
3090 y1y=ry(jp,i)/dis(i,jp)
3091 y1z=rz(jp,i)/dis(i,jp)
3096 yy33=1.0D0/(dis(i,jp)*dis(i,ip))
3097 yyy3a=pin1(im,j)/(dis(i,jp)**2)
3098 yyy3b=pin1(im,j)/(dis(i,ip)**2)
3099 yy3=-pin1(im,j)/dshe
3100 y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
3101 $ +yyy3b*rx(i,ip))*yy3
3102 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
3103 $ +yyy3b*ry(i,ip))*yy3
3104 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
3105 $ +yyy3b*rz(i,ip))*yy3
3107 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3108 yyy4=pin2(im,j)/(dis(i,jp)**2)
3109 yy4=-pin2(im,j)/dshe
3110 y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3111 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3112 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3114 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3115 yyy5=pin3(im,j)/(dis(i,ip)**2)
3116 yy5=-pin3(im,j)/dshe
3117 y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3118 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3119 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3132 shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3133 $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3134 shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3135 $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3136 shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3137 $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3138 ! shefx(i,6)=shefx(i,6)
3139 ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3140 ! shefy(i,6)=shefy(i,6)
3141 ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3142 ! shefz(i,6)=shefz(i,6)
3143 ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3145 yy6=-(dis(jpp,i)-uldhb)/dldhb
3146 y6x=rx(jpp,i)/dis(jpp,i)
3147 y6y=ry(jpp,i)/dis(jpp,i)
3148 y6z=rz(jpp,i)/dis(jpp,i)
3153 yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3154 yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3155 yyy8b=pina1(im,j)/(dis(i,ip)**2)
3156 yy8=-pina1(im,j)/dshe
3157 y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3158 $ +yyy8b*rx(i,ip))*yy8
3159 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3160 $ +yyy8b*ry(i,ip))*yy8
3161 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3162 $ +yyy8b*rz(i,ip))*yy8
3164 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3165 yyy9=pina2(im,j)/(dis(i,jpp)**2)
3166 yy9=-pina2(im,j)/dshe
3167 y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3168 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3169 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3171 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3172 yyy10=pina3(im,j)/(dis(i,ip)**2)
3173 yy10=-pina3(im,j)/dshe
3174 y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3175 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3176 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3178 sx=y66x+y8x+y9x+y10x
3179 sy=y66y+y8y+y9y+y10y
3180 sz=y66z+y8z+y9z+y10z
3189 shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3190 $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3191 shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3192 $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3193 shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3194 $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3196 ! shefx(i,6)=shefx(i,6)
3197 ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3198 ! shefy(i,6)=shefy(i,6)
3199 ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3200 ! shefz(i,6)=shefz(i,6)
3201 ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3211 c-----------------------------------------------------------------------
3212 subroutine sheetforce11
3215 parameter(maxca=800)
3216 real*8 dfa_cutoff,dfa_cutoff_delta
3217 parameter(dfa_cutoff=15.5d0)
3218 parameter(dfa_cutoff_delta=0.5d0)
3219 cc**********************************************************************
3220 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3221 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3222 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3223 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3224 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3225 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3226 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3227 real*8 rx(maxca,maxca)
3228 real*8 ry(maxca,maxca),rz(maxca,maxca)
3229 real*8 bx(maxca),by(maxca),bz(maxca)
3230 real*8 dis(maxca,maxca)
3231 real*8 shefx(maxca,12),shefy(maxca,12)
3232 real*8 shefz(maxca,12)
3233 real*8 dp45,dm45,w_beta
3234 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3235 $ c00,s00,ulnex,dnex
3236 integer inb,nmax,iselect
3237 cc**********************************************************************
3238 common /phys1/ inb,nmax,iselect
3240 common /difvec/ rx,ry,rz
3241 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3242 $ c00,s00,ulnex,dnex
3243 common /sheconst/ dp45,dm45,w_beta
3244 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3245 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3246 common /shef/ shefx,shefy,shefz
3247 ci integer istrand(maxca,maxca)
3248 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3249 ci common /shetest/ istrand,istrand_p,istrand_m
3250 C********************************************************************************
3252 integer j,jm,jmm,ip,i,ipp
3253 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3254 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3255 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3256 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3257 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3258 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3259 real*8 e_gcont,fprim_gcont
3260 C********************************************************************************
3267 if (dis(i,jmm).lt.dfa_cutoff) then
3268 call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3269 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3274 ci if(istrand(i,jmm).eq.1
3275 ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3278 yy1=-(dis(ipp,j)-ulhb)/dlhb
3279 y1x=rx(ipp,j)/dis(ipp,j)
3280 y1y=ry(ipp,j)/dis(ipp,j)
3281 y1z=rz(ipp,j)/dis(ipp,j)
3286 yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3287 yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3288 yy3=-pin2(i,jmm)/dshe
3289 y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3290 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3291 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3293 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3294 yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3295 yy4=-pin3(i,jmm)/dshe
3296 y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3297 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3298 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3300 yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3301 yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3302 yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3303 yy5=-pin4(i,jmm)/dshe
3304 y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3305 $ -yyy5b*rx(jm,j))*yy5
3306 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3307 $ -yyy5b*ry(jm,j))*yy5
3308 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3309 $ -yyy5b*rz(jm,j))*yy5
3322 shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3323 $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3324 shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3325 $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3326 shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3327 $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3329 ! shefx(j,11)=shefx(j,11)
3330 ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3331 ! shefy(j,11)=shefy(j,11)
3332 ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3333 ! shefz(j,11)=shefz(j,11)
3334 ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3336 yy6=-(dis(ip,j)-uldhb)/dldhb
3337 y6x=rx(ip,j)/dis(ip,j)
3338 y6y=ry(ip,j)/dis(ip,j)
3339 y6z=rz(ip,j)/dis(ip,j)
3344 yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3345 yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3346 yy8=-pina1(i,jmm)/dshe
3347 y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3348 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3349 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3351 yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3352 yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3353 yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3354 yy9=-pina2(i,jmm)/dshe
3355 y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3356 $ -yyy9b*rx(jm,j))*yy9
3357 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3358 $ -yyy9b*ry(jm,j))*yy9
3359 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3360 $ -yyy9b*rz(jm,j))*yy9
3362 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3363 yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3364 yy10=-pina4(i,jmm)/dshe
3365 y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3366 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3367 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3369 sx=y66x+y8x+y9x+y10x
3370 sy=y66y+y8y+y9y+y10y
3371 sz=y66z+y8z+y9z+y10z
3380 shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3381 $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3382 shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3383 $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3384 shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3385 $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3387 ! shefx(j,11)=shefx(j,11)
3388 ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3389 ! shefy(j,11)=shefy(j,11)
3390 ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3391 ! shefz(j,11)=shefz(j,11)
3392 ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3402 c-----------------------------------------------------------------------
3403 subroutine sheetforce12
3406 parameter(maxca=800)
3407 real*8 dfa_cutoff,dfa_cutoff_delta
3408 parameter(dfa_cutoff=15.5d0)
3409 parameter(dfa_cutoff_delta=0.5d0)
3410 cc**********************************************************************
3411 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3412 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3413 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3414 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3415 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3416 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3417 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3418 real*8 rx(maxca,maxca)
3419 real*8 ry(maxca,maxca),rz(maxca,maxca)
3420 real*8 bx(maxca),by(maxca),bz(maxca)
3421 real*8 dis(maxca,maxca)
3422 real*8 shefx(maxca,12),shefy(maxca,12)
3423 real*8 shefz(maxca,12)
3424 real*8 dp45,dm45,w_beta
3425 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3426 $ c00,s00,ulnex,dnex
3427 integer inb,nmax,iselect
3428 cc**********************************************************************
3429 common /phys1/ inb,nmax,iselect
3431 common /difvec/ rx,ry,rz
3432 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3433 $ c00,s00,ulnex,dnex
3434 common /sheconst/ dp45,dm45,w_beta
3435 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3436 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3437 common /shef/ shefx,shefy,shefz
3438 ci integer istrand(maxca,maxca)
3439 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3440 ci common /shetest/ istrand,istrand_p,istrand_m
3441 cc**********************************************************************
3443 integer j,jm,jmm,ip,i,ipp,jp
3444 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3445 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3446 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3447 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3448 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3449 real*8 e_gcont,fprim_gcont
3450 !c*************************************************************************c
3456 if (dis(i,jm).lt.dfa_cutoff) then
3457 call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3458 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3463 ci if(istrand(i,jm).eq.1
3464 ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3467 yy1=-(dis(ip,j)-ulhb)/dlhb
3468 y1x=rx(ip,j)/dis(ip,j)
3469 y1y=ry(ip,j)/dis(ip,j)
3470 y1z=rz(ip,j)/dis(ip,j)
3475 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3476 yyy3=pin1(i,jm)/(dis(ip,j)**2)
3477 yy3=-pin1(i,jm)/dshe
3478 y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3479 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3480 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3481 yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3483 yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3484 yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3485 yy4=-pin2(i,jm)/dshe
3486 y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3487 $ +yyy4b*rx(j,jp))*yy4
3488 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3489 $ +yyy4b*ry(j,jp))*yy4
3490 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3491 $ +yyy4b*rz(j,jp))*yy4
3493 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3494 yyy5=pin4(i,jm)/(dis(j,jp)**2)
3495 yy5=-pin4(i,jm)/dshe
3496 y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3497 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3498 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3511 shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3512 $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3513 shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3514 $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3515 shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3516 $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3518 ! shefx(j,12)=shefx(j,12)
3519 ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3520 ! shefy(j,12)=shefy(j,12)
3521 ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3522 ! shefz(j,12)=shefz(j,12)
3523 ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3525 yy6=-(dis(ipp,j)-uldhb)/dldhb
3526 y6x=rx(ipp,j)/dis(ipp,j)
3527 y6y=ry(ipp,j)/dis(ipp,j)
3528 y6z=rz(ipp,j)/dis(ipp,j)
3533 yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3534 yyy8=pina2(i,jm)/(dis(j,jp)**2)
3535 yy8=-pina2(i,jm)/dshe
3536 y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3537 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3538 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3540 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3541 yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3542 yy9=-pina3(i,jm)/dshe
3543 y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3544 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3545 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3547 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3548 yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3549 yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3550 yy10=-pina4(i,jm)/dshe
3551 y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3552 $ +yyy10b*rx(j,jp))*yy10
3553 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3554 $ +yyy10b*ry(j,jp))*yy10
3555 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3556 $ +yyy10b*rz(j,jp))*yy10
3558 sx=y66x+y8x+y9x+y10x
3559 sy=y66y+y8y+y9y+y10y
3560 sz=y66z+y8z+y9z+y10z
3569 shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3570 $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3571 shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3572 $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3573 shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3574 $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3585 C===============================================================================