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 'DIMENSIONS.ZSCOPT'
75 include 'COMMON.IOUNITS'
76 include 'COMMON.CHAIN'
78 include 'COMMON.FFIELD'
81 C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
82 C THIS SHOULD BE MODIFIED!!
89 integer ica1, ica2,ica3,ica4,ica5
90 integer ishell, inca, itmp,iitmp
95 open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
97 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
100 write(iout,'(a)') 'dist_dfa.dat is opened!'
102 read(iodfa, '(a)') buffer
103 C read number of restraints
104 read(iodfa, *) IDFADIS
105 read(iodfa, *) dis_inc
107 read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
126 C READ ANGLE RESTRAINTS
128 open(iodfa, file='phi_dfa.dat',status='old',err=35)
130 35 write(iout,'(a)') 'Error opening dist_dfa.dat file'
134 write(iout,'(a)') 'phi_dfa.dat is opened!'
137 read(iodfa, '(a)') buffer
138 C READ NUMBER OF RESTRAINTS
139 READ(iodfa, *) IDFAPHI
140 read(iodfa,*) phi_inc
142 read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
153 read(iodfa,*) tmp1,tmp2
167 open(iodfa, file='theta_dfa.dat',status='old',err=41)
169 41 write(iout,'(a)') 'Error opening dist_dfa.dat file'
172 write(iout,'(a)') 'theta_dfa.dat is opened!'
174 read(iodfa, '(a)') buffer
175 C READ NUMBER OF RESTRAINTS
176 READ(iodfa, *) IDFATHE
177 read(iodfa,*) the_inc
180 read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
191 read(iodfa,*) tmp1,tmp2
203 C END of READING ANGLE RESTRAINT!
205 C NUMBER OF NEIGHBOR CAs
206 open(iodfa,file='nei_dfa.dat',status='old',err=37)
208 37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
211 write(iout,'(a)') 'nei_dfa.dat is opened!'
213 read(iodfa, '(a)') buffer
214 C READ NUMBER OF RESTRAINTS
215 READ(iodfa, *) idfanei
216 read(iodfa,*) nei_inc
219 read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
228 C write(*,*) 'READ NEI:',i,j,fnei(i,j)
238 C END OF NEIGHBORING CA
240 C READ BETA RESTRAINT
241 if (wdfa_beta.eq.0.0) return
242 open(iodfa, file='beta_dfa.dat',status='old',err=39)
244 39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
247 write(iout,'(a)') 'beta_dfa.dat is opened!'
249 read(iodfa,'(a)') buffer
251 read(iodfa,*) beta_inc
254 read(iodfa,*) ica1, iitmp
258 c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
263 C END OF BETA RESTRAINT
268 subroutine edfad(edfadis)
270 implicit real*8 (a-h,o-z)
272 include 'COMMON.CHAIN'
273 include 'COMMON.DERIV'
276 double precision edfadis
277 integer i, iatm1, iatm2,idiff
278 double precision ckk, sckk,dist,texp
279 double precision jix,jiy,jiz,ep,fp,scc
286 iatm1=idislis(1,i)+ishiftca
287 iatm2=idislis(2,i)+ishiftca
288 idiff = abs(iatm1-iatm2)
290 JIX=c(1,iatm2)-c(1,iatm1)
291 JIY=c(2,iatm2)-c(2,iatm1)
292 JIZ=c(3,iatm2)-c(3,iatm1)
293 DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
306 if (dtmp.ge.15.0d0) then
309 c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
310 texp = exp(-dtmp)/sckk
313 ep=ep+sccdist(i,j)*texp
314 fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
316 C write(*,'(2i8,6f12.5)') i, j, dist,
317 C & fdist(i,j), ep, fp, sccdist(i,j), scc
325 c IF(ABS(EP).lt.1.0d-20)THEN
328 c IF (ABS(FP).lt.1.0d-20) THEN
332 edfadis=edfadis+ep*dis_inc*wwdist
334 gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
335 gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
336 gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
338 gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
339 gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
340 gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
347 subroutine edfat(edfator)
349 implicit real*8 (a-h,o-z)
351 include 'COMMON.CHAIN'
352 include 'COMMON.DERIV'
357 double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
358 double precision cwidth, cwidth2
359 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
360 parameter (TENM20=1.0d-20)
372 iatom(iii)=iphilis(iii,i)+ishiftca
375 C ANGLE VECTOR CALCULTION
376 RIX=C(1,IATOM(2))-C(1,IATOM(1))
377 RIY=C(2,IATOM(2))-C(2,IATOM(1))
378 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
380 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
381 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
382 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
384 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
385 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
386 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
388 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
389 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
390 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
392 GIX=RIY*RIPZ-RIZ*RIPY
393 GIY=RIZ*RIPX-RIX*RIPZ
394 GIZ=RIX*RIPY-RIY*RIPX
396 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
397 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
398 GIPZ=RIPX*RIPPY-RIPY*RIPPX
400 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
401 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
402 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
404 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
405 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
406 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
408 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
409 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
410 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
412 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
413 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
414 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
415 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
417 C END OF ANGLE VECTOR CALCULTION
419 TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
420 APHI(1)=TDOT/(DGI*DRIPP)
421 TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
422 APHI(2)=TDOT/(DGIP*DRIP3)
430 DDPS1=APHI(1)-FPHI1(i,j)
431 DDPS2=APHI(2)-FPHI2(i,j)
433 DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
435 if (dtmp.ge.15.0d0) then
438 c ps_tmp = dfaexp(idint(dtmp*1000)+1)
442 ephi=ephi+sccphi(i,j)*ps_tmp
444 tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
445 tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
448 c write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
449 c & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
452 ephi=-ephi/scc*phi_inc*wwangle
453 tfphi1=tfphi1/scc*phi_inc*wwangle
454 tfphi2=tfphi2/scc*phi_inc*wwangle
456 IF (ABS(EPHI).LT.1d-20) THEN
459 IF (ABS(TFPHI1).LT.1d-20) THEN
462 IF (ABS(TFPHI2).LT.1d-20) THEN
466 C FORCE DIRECTION CALCULATION
471 DM1=1.0d0/(DGI*DRIPP)
473 GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
474 DM2=GIRPP/(DGI**3*DRIPP)
475 DM3=GIRPP/(DGI*DRIPP**3)
477 DM4=1.0d0/(DGIP*DRIP3)
479 GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
480 DM5=GIRP3/(DGIP**3*DRIP3)
481 DM6=GIRP3/(DGIP*DRIP3**3)
483 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
484 & +( GIZ* RIPY- GIY* RIPZ)*DM2
485 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
486 & +( GIX* RIPZ- GIZ* RIPX)*DM2
487 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
488 & +( GIY* RIPX- GIX* RIPY)*DM2
492 C SECOND ATOM BY PHI1
493 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
494 & -(CIPY*GIZ-CIPZ*GIY)*DM2
495 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
496 & -(CIPZ*GIX-CIPX*GIZ)*DM2
497 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
498 & -(CIPX*GIY-CIPY*GIX)*DM2
502 C SECOND ATOM BY PHI2
504 & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
505 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
507 & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
508 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
510 & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
511 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
513 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
514 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
515 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
516 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
517 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
518 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
524 & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
525 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
527 & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
528 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
530 & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
531 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
532 C FOURTH ATOM BY PHI1
533 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
534 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
535 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
536 C FOURTH ATOM BY PHI2
538 & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
539 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
540 & + RIP3X*DM6)*TFPHI2
542 & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
543 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
544 & + RIP3Y*DM6)*TFPHI2
546 & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
547 & -( GIPX*RIPY-RIPX*GIPY)*DM5
548 & + RIP3Z*DM6)*TFPHI2
550 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
551 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
552 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
553 C END OF FORCE DIRECTION
556 gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
557 gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
558 gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
561 enephi = enephi + ephi
562 c end of single assignment statement
564 C END OF PHI RESTRAINT
566 C START OF THETA ANGLE
571 iatom(iii)=ithelis(iii,i)+ishiftca
575 C ANGLE VECTOR CALCULTION
576 RIX=C(1,IATOM(2))-C(1,IATOM(1))
577 RIY=C(2,IATOM(2))-C(2,IATOM(1))
578 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
580 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
581 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
582 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
584 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
585 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
586 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
588 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
589 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
590 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
592 GIX=RIY*RIPZ-RIZ*RIPY
593 GIY=RIZ*RIPX-RIX*RIPZ
594 GIZ=RIX*RIPY-RIY*RIPX
596 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
597 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
598 GIPZ=RIPX*RIPPY-RIPY*RIPPX
600 GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
601 GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
602 GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
604 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
605 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
606 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
608 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
609 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
610 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
612 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
613 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
614 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
616 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
617 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
618 DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
619 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
620 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
621 C END OF ANGLE VECTOR CALCULTION
623 TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
624 ATHE(1)=TDOT/(DGI*DGIP)
625 TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
626 ATHE(2)=TDOT/(DGIP*DGIPP)
635 ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
636 ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
637 dtmp= (ddth1**2+ddth2**2)/cwidth2
638 if ( dtmp .ge. 15.0d0) then
641 c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
645 ethe=ethe+sccthe(i,j)*th_tmp
647 tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
648 tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
650 c write(2,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
651 c & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
654 ethe=-ethe/scc*the_inc*wwangle
655 tfthe1=tfthe1/scc*the_inc*wwangle
656 tfthe2=tfthe2/scc*the_inc*wwangle
658 IF (ABS(ETHE).LT.TENM20) THEN
661 IF (ABS(TFTHE1).LT.TENM20) THEN
664 IF (ABS(TFTHE2).LT.TENM20) THEN
673 DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
674 DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
676 DM4=1.0d0/(DGIP*DGIPP)
677 DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
678 DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
680 C FIRST ATOM BY THETA1
681 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
682 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
683 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
684 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
685 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
686 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
687 C SECOND ATOM BY THETA1
688 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
689 & -(CIPY*GIZ-CIPZ*GIY)*DM2
690 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
691 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
692 & -(CIPZ*GIX-CIPX*GIZ)*DM2
693 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
694 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
695 & -(CIPX*GIY-CIPY*GIX)*DM2
696 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
697 C SECOND ATOM BY THETA2
699 & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
700 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
702 & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
703 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
705 & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
706 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
707 C THIRD ATOM BY THETA1
708 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
709 & -(GIY*RIZ-GIZ*RIY)*DM2
710 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
711 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
712 & -(GIZ*RIX-GIX*RIZ)*DM2
713 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
714 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
715 & -(GIX*RIY-GIY*RIX)*DM2
716 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
717 C THIRD ATOM BY THETA2
719 & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
720 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
721 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
723 & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
724 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
725 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
727 & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
728 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
729 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
730 C FOURTH ATOM BY THETA1
731 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
732 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
733 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
734 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
735 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
736 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
737 C FOURTH ATOM BY THETA2
739 & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
740 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
741 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
743 & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
744 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
745 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
747 & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
748 & -(GIPX*RIPY-GIPY*RIPX)*DM5
749 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
750 C FIFTH ATOM BY THETA2
751 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
752 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
753 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
754 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
755 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
756 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
757 C !! END OF FORCE DIRECTION!!!!
759 gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
760 gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
761 gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
764 enethe = enethe + ethe
767 edfator = enephi + enethe
772 subroutine edfan(edfanei)
773 C DFA neighboring CA restraint
774 implicit real*8 (a-h,o-z)
776 include 'COMMON.CHAIN'
777 include 'COMMON.DERIV'
781 integer kshnum, n1atom
783 double precision enenei,tmp_n
784 double precision pai,hpai
785 double precision jix,jiy,jiz,ndiff,snorm_nei
786 double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
787 double precision dr,dr2,half,ntmp,dtmp
789 parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
790 parameter(pai=3.14159265358979323846D0)
791 parameter(hpai=1.5707963267948966D0)
792 parameter(snorm_nei=0.886226925452758D0)
798 c print*, 's1:', s1(:)
799 c print*, 's2:', s2(:)
804 n1atom=ineilis(i)+ishiftca
805 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
818 do j = ishiftca+1, ilastca
820 if (n1atom.eq.j) cycle
822 jix=c(1,j)-c(1,n1atom)
823 jiy=c(2,j)-c(2,n1atom)
824 jiz=c(3,j)-c(3,n1atom)
825 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
827 c write(*,*) n1atom, j, dist
830 if (dist.lt.s1(kshnum).and.
831 & dist.gt.s2(kshnum-1)) then
835 c write(*,*) 'case1:',tmp_n
844 elseif(dist.ge.s1(kshnum).and.
845 & dist.le.s2(kshnum)) then
847 dnei=(dist-s1(kshnum))/dr2*pai
848 tmp_n=tmp_n + half*(1+cos(dnei))
849 c write(*,*) 'case2:',tmp_n
850 ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
860 elseif(dist.ge.s1(kshnum-1).and.
861 & dist.le.s2(kshnum-1)) then
862 dnei=(dist-s1(kshnum-1))/dr2*pai
863 tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
864 c write(*,*) 'case3:',tmp_n
865 ftmp = hpai*sin(dnei)/dr2/dist
877 elseif(kshnum.eq.1) then
879 if(dist.lt.s1(kshnum))then
882 c write(*,*) 'case4:',tmp_n
890 elseif(dist.ge.s1(kshnum).and.
891 & dist.le.s2(kshnum))then
893 dnei=(dist-s1(kshnum))/dr2*pai
894 tmp_n=tmp_n + half*(1+cos(dnei))
895 c write(*,*) 'case5:',tmp_n
896 ftmp = -hpai*sin(dnei)/dr2/dist
917 ndiff = tmp_n-fnei(i,imin)
920 if (dtmp.ge.15.0d0) then
923 c ntmp = dfaexp( idint(dtmp*1000) + 1 )
927 enei=enei+sccnei(i,imin)*ntmp
929 & sccnei(i,imin)*ntmp*ndiff*2.0d0
930 scc=scc+sccnei(i,imin)
932 c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
933 c & fnei(i,imin),sccnei(i,imin),enei,scc
936 enei=-enei/scc*snorm_nei*nei_inc*wwnei
937 tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
939 c if (abs(enei).lt.1.0d-20)then
942 c if (abs(tmp_fnei).lt.1.0d-20) then
951 do j=ishiftca+1,ilastca
952 t2dx(j)=t2dx(j)*tmp_fnei
953 t2dy(j)=t2dy(j)*tmp_fnei
954 t2dz(j)=t2dz(j)*tmp_fnei
957 gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
958 gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
959 gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
961 do j=ishiftca+1,ilastca
962 gdfan(1,j)=gdfan(1,j)+t2dx(j)
963 gdfan(2,j)=gdfan(2,j)+t2dy(j)
964 gdfan(3,j)=gdfan(3,j)+t2dz(j)
977 subroutine edfab(edfabeta)
979 implicit real*8 (a-h,o-z)
982 include 'COMMON.CHAIN'
983 include 'COMMON.DERIV'
987 parameter(PAI=3.14159265358979323846D0)
988 parameter (maxca=800)
990 real*8 bx(maxres),by(maxres),bz(maxres)
991 real*8 vbet(maxres,maxres)
992 real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
993 real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
994 real*8 vbeta,vbetp,vbetm
995 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
997 real*8 dp45,dm45,w_beta
999 real*8 cph(maxca),cth(maxca)
1000 real*8 atx(maxca),aty(maxca),atz(maxca)
1001 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1002 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1003 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1005 real*8 astx(maxca),asty(maxca),astz(maxca)
1006 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1007 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1008 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1010 real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
1011 & astxnum(maxca),astynum(maxca),astznum(maxca),
1012 & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
1013 & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
1014 & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
1015 & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
1016 & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
1017 & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
1018 & cth_orig(maxca),sth_orig(maxca)
1020 common /sheca/ bx,by,bz
1021 common /shee/ vbeta,vbet,vbetp,vbetm
1022 common /shetf/ shetfx,shetfy,shetfz
1023 common /shef/ shefx, shefy, shefz
1024 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
1025 & c00,s00,ulnex,dnex
1026 common /sheconst/ dp45,dm45,w_beta
1028 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1029 $ atmmz,atm3x,atm3y,atm3z
1030 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1031 $ astmmz,astm3x,astm3y,astm3z
1033 common /coscos/ cph,cth
1036 C End of sheet variables
1039 double precision enebet
1042 c bx=0.0d0;by=0.0d0;bz=0.0d0
1043 c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1047 do i=ishiftca+1,ilastca
1048 bx(i-ishiftca)=c(1,i)
1049 by(i-ishiftca)=c(2,i)
1050 bz(i-ishiftca)=c(3,i)
1053 c do i=1,ilastca-ishiftca
1054 c read(99,*) bx(i),by(i),bz(i)
1062 ULNEX=COS(60.0D0/180.0D0*PAI)
1069 C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1070 S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1076 C END OF INITIALIZATION
1078 nca=ilastca-ishiftca
1080 call angvectors(nca)
1081 call sheetforce(nca,wshet)
1083 c end of sheet energy and force
1086 shetfx(j)=shetfx(j)*beta_inc
1087 shetfy(j)=shetfy(j)*beta_inc
1088 shetfz(j)=shetfz(j)*beta_inc
1089 c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1092 vbeta=vbeta*beta_inc
1097 gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1098 gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1099 gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1104 write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1118 call angvectors(nca)
1120 call angvectors(nca)
1121 atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1122 astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1124 atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1125 astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1128 atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1129 astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1132 atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1133 astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1137 call angvectors(nca)
1139 call angvectors(nca)
1141 atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1142 astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1144 atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1145 astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1148 atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1149 astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1152 atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1153 astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1157 call angvectors(nca)
1159 call angvectors(nca)
1162 atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1163 astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1165 atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1166 astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1169 atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1170 astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1173 atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1174 astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1180 write (*,'(2i5,a2,6f10.5)')
1181 & i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1182 & astxnum(i),astx(i),astxnum(i)/astx(i),
1183 & i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1184 & astynum(i),asty(i),astynum(i)/asty(i),
1185 & i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1186 & astznum(i),astz(i),astznum(i)/astz(i),
1187 & i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1188 & astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1189 & i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1190 & astmynum(i),astmy(i),astmynum(i)/astmy(i),
1191 & i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1192 & astmznum(i),astmz(i),astmznum(i)/astmz(i),
1193 & i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1194 & astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1195 & i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1196 & astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1197 & i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1198 & astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1199 & i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1200 & astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1201 & i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1202 & astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1203 & i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1204 & astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1205 & i,0," ",cth_orig(i),sth_orig(i)
1215 call angvectors(nca)
1216 call sheetforce(nca,wshet)
1217 vbeta1=vbeta*beta_inc
1219 call angvectors(nca)
1220 call sheetforce(nca,wshet)
1221 vbeta2=vbeta*beta_inc
1222 gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1226 call angvectors(nca)
1227 call sheetforce(nca,wshet)
1228 vbeta1=vbeta*beta_inc
1230 call angvectors(nca)
1231 call sheetforce(nca,wshet)
1232 vbeta2=vbeta*beta_inc
1233 gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1237 call angvectors(nca)
1238 call sheetforce(nca,wshet)
1239 vbeta1=vbeta*beta_inc
1241 call angvectors(nca)
1242 call sheetforce(nca,wshet)
1243 vbeta2=vbeta*beta_inc
1244 gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1251 call angvectors(nca)
1252 call sheetforce(nca,wshet)
1254 shetfx(j)=shetfx(j)*beta_inc
1255 shetfy(j)=shetfy(j)*beta_inc
1256 shetfz(j)=shetfz(j)*beta_inc
1260 write(*,*) 'xyz analytical and numerical gradient'
1262 write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1263 & ,(-gdfab(i,j),i=1,3)
1267 write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1268 & shetfy(j)/gdfab(2,j),
1269 & shetfz(j)/gdfab(3,j)
1277 C-------------------------------------------------------------------------------
1278 subroutine angvectors(nca)
1279 c implicit real*4(a-h,o-z)
1283 parameter(maxca=800)
1285 parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1287 real*8 bx(maxca),by(maxca),bz(maxca)
1288 real*8 dis(maxca,maxca)
1289 real*8 apx(maxca),apy(maxca),apz(maxca)
1290 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
1291 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
1292 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
1293 real*8 atx(maxca),aty(maxca),atz(maxca)
1294 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1295 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1296 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1297 real*8 astx(maxca),asty(maxca),astz(maxca)
1298 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1299 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1300 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1302 real*8 cph(maxca),cth(maxca)
1305 integer i, ip, ipp, ip3, j
1306 real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1307 real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1308 real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1309 real*8 cix, ciy, ciz, cipx, cipy, cipz
1310 real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1311 real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
1312 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1313 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1314 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1315 real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1316 real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1317 real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1318 real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1319 real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1320 real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1321 real*8 grpp,gx,gy,gz
1322 real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1323 real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1324 integer inb,nmax,iselect
1326 common /sheca/ bx,by,bz
1327 common /difvec/ rx, ry, rz
1328 common /ulang/ ulcos
1329 common /phys1/ inb,nmax,iselect
1332 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1333 & apmmz,apm3x,apm3y,apm3z
1334 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1335 & atmmz,atm3x,atm3y,atm3z
1336 common /coscos/ cph,cth
1337 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1338 & astmmz,astm3x,astm3y,astm3z
1340 C-------------------------------------------------------------------------------
1341 c write(*,*) 'inside angvectors'
1346 cph=zero; cth=zero; sth=zero
1347 apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1348 apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1349 atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1350 atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1351 astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1352 astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1355 C r[x,y,z] calc and distance calculation
1356 rx=zero;ry=zero;rz=zero
1363 dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1364 c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1365 c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1366 c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1367 c write(*,*) 'dis(i,j):',i,j,dis(i,j)
1370 c end of r[x,y,z] calc
1376 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1377 ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1378 $ +rz(i,ip)*rz(ip,ipp)
1379 ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1382 c end of virtual bond angle
1383 c write(*,*) 'inside angvectors1'
1395 rippx=bx(ip3)-bx(ipp)
1396 rippy=by(ip3)-by(ipp)
1397 rippz=bz(ip3)-bz(ipp)
1399 gx=riy*ripz-riz*ripy
1400 gy=riz*ripx-rix*ripz
1401 gz=rix*ripy-riy*ripx
1402 gpx=ripy*rippz-ripz*rippy
1403 gpy=ripz*rippx-ripx*rippz
1404 gpz=ripx*rippy-ripy*rippx
1405 gpcrp_x=gpy*ripz-gpz*ripy
1406 gpcrp_y=gpz*ripx-gpx*ripz
1407 gpcrp_z=gpx*ripy-gpy*ripx
1408 d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1409 gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1410 & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1416 gmx=rimy*riz-rimz*riy
1417 gmy=rimz*rix-rimx*riz
1418 gmz=rimx*riy-rimy*rix
1419 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1421 ggm=gmx*gx+gmy*gy+gmz*gz
1422 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1428 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1430 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1431 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1433 c write(*,*) 'inside angvectors2'
1435 rimmx=bx(i-1)-bx(i-2)
1436 rimmy=by(i-1)-by(i-2)
1437 rimmz=bz(i-1)-bz(i-2)
1439 gmmx=rimmy*rimz-rimmz*rimy
1440 gmmy=rimmz*rimx-rimmx*rimz
1441 gmmz=rimmx*rimy-rimmy*rimx
1442 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1444 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1445 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1446 gmcrim_x=gmy*rimz-gmz*rimy
1447 gmcrim_y=gmz*rimx-gmx*rimz
1448 gmcrim_z=gmx*rimy-gmy*rimx
1449 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1450 d_gmcrim3=d_gmcrim**3
1451 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1452 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1456 rim3x=bx(i-2)-bx(i-3)
1457 rim3y=by(i-2)-by(i-3)
1458 rim3z=bz(i-2)-bz(i-3)
1459 g3x=rim3y*rimmz-rim3z*rimmy
1460 g3y=rim3z*rimmx-rim3x*rimmz
1461 g3z=rim3x*rimmy-rim3y*rimmx
1462 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1463 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1464 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1465 cc**********************************************************************
1466 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1467 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1468 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1469 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1470 d_gmmcrimm3=d_gmmcrimm**3
1471 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1472 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1479 dg=sqrt(gx**2+gy**2+gz**2)
1480 dgp=sqrt(gpx**2+gpy**2+gpz**2)
1483 ggp=gx*gpx+gy*gpy+gz*gpz
1484 grpp=gx*rippx+gy*rippy+gz*rippz
1486 if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1487 & .and.d_gpcrp.gt.0.0D0) then
1488 cph(i)=grpp/dg/dripp
1490 sth(i)=gpcrp__g/d_gpcrp/dg
1498 c write(*,*) 'inside angvectors3'
1500 if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1501 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1504 d12=1.0D0/(dg*dripp)
1505 d13=grpp/(dg3*dripp)
1506 sd10=1.0D0/(d_gpcrp*dg)
1507 sd11=gpcrp__g/(d_gpcrp*dg3)
1517 atx(i)=(ripz*gpy-ripy*gpz)*d10
1518 & -(gy*ripz-gz*ripy)*d11
1519 aty(i)=(ripx*gpz-ripz*gpx)*d10
1520 & -(gz*ripx-gx*ripz)*d11
1521 atz(i)=(ripy*gpx-ripx*gpy)*d10
1522 & -(gx*ripy-gy*ripx)*d11
1523 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1524 & +ripy*gpy*ripx-gpx*ripz**2)
1525 & -sd11*(gy*ripz-gz*ripy)
1526 asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1527 & -gpy*ripx**2+gpz*ripy*ripz)
1528 & -sd11*(-gx*ripz+gz*ripx)
1529 astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1530 & -gpz*ripy**2+ripz*gpx*ripx)
1531 & -sd11*(gx*ripy-gy*ripx)
1532 apx(i)=(ripz*rippy-ripy*rippz)*d12
1533 & -(gy*ripz-gz*ripy)*d13
1534 apy(i)=(ripx*rippz-ripz*rippx)*d12
1535 & -(gz*ripx-gx*ripz)*d13
1536 apz(i)=(ripy*rippx-ripx*rippy)*d12
1537 & -(gx*ripy-gy*ripx)*d13
1549 if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1550 & .and.d_gcr3.gt.0.0D0) then
1554 d23=1.0D0/(dgm*drip)
1555 d24=gmrp/(dgm3*drip)
1556 sd20=1.0D0/(d_gcr*dgm)
1557 sd21=gcr__gm/(d_gcr3*dgm)
1558 sd22=gcr__gm/(d_gcr*dgm3)
1569 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1570 & -(ciy*gmz-ciz*gmy)*d21
1571 & +(ripy*gz-ripz*gy)*d22
1572 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1573 & -(ciz*gmx-cix*gmz)*d21
1574 & +(ripz*gx-ripx*gz)*d22
1575 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1576 & -(cix*gmy-ciy*gmx)*d21
1577 & +(ripx*gy-ripy*gx)*d22
1578 cc**********************************************************************
1579 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1580 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1581 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1582 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1583 & +gcr_z*(-ripz*rix+gy))
1584 & -sd22*(-gmy*ciz+gmz*ciy)
1586 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1587 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1589 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1590 & -gcr_z*(ripz*riy+gx))
1591 & -sd22*(gmx*ciz-gmz*cix)
1593 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1594 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1596 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1597 & +gcr_z*(ripy*riy+ripx*rix))
1598 & -sd22*(-gmx*ciy+gmy*cix)
1599 cc**********************************************************************
1600 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1601 & -(ciy*gmz-ciz*gmy)*d24
1602 apmy(i)=(ciz*ripx-ripz*cix)*d23
1603 & -(ciz*gmx-cix*gmz)*d24
1604 apmz(i)=(cix*ripy-ripx*ciy)*d23
1605 & -(cix*gmy-ciy*gmx)*d24
1609 if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1610 & .and.d_gmcrim3.gt.0.0D0) then
1611 d30=1.0D0/(dgm*dgmm)
1612 d31=gmmgm/(dgm3*dgmm)
1613 d32=gmmgm/(dgm*dgmm3)
1614 d33=1.0D0/(dgmm*dri)
1615 d34=gmmr/(dgmm3*dri)
1616 d35=gmmr/(dgmm*dri3)
1617 sd30=1.0D0/(d_gmcrim*dgmm)
1618 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1619 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1632 c write(*,*) 'inside angvectors4'
1634 cc**********************************************************************
1635 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1636 & -(ciy*gmz-ciz*gmy)*d31
1637 & -(gmmy*rimmz-gmmz*rimmy)*d32
1638 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1639 & -(ciz*gmx-cix*gmz)*d31
1640 & -(gmmz*rimmx-gmmx*rimmz)*d32
1641 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1642 & -(cix*gmy-ciy*gmx)*d31
1643 & -(gmmx*rimmy-gmmy*rimmx)*d32
1644 cc**********************************************************************
1645 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1646 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1647 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1648 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1649 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1650 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1651 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1653 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1654 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1655 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1656 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1657 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1658 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1659 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1661 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1662 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1663 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1664 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1665 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1666 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1667 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1668 c**********************************************************************
1669 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1670 & -(gmmy*rimmz-gmmz*rimmy)*d34
1672 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1673 & -(gmmz*rimmx-gmmx*rimmz)*d34
1675 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1676 & -(gmmx*rimmy-gmmy*rimmx)*d34
1681 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1682 & .and.drim3.gt.0.0D0
1683 & .and.d_gmmcrimm3.gt.0.0D0) then
1684 d40=1.0D0/(dg30*dgmm)
1685 d41=g3gmm/(dg30*dgmm3)
1686 d42=1.0D0/(dg30*drim)
1687 d43=g3rim_/(dg30*drim3)
1688 sd40=1.0D0/(dg30*d_gmmcrimm)
1689 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1698 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1699 & -(gmmy*rimmz-gmmz*rimmy)*d41
1700 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1701 & -(gmmz*rimmx-gmmx*rimmz)*d41
1702 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1703 & -(gmmx*rimmy-gmmy*rimmx)*d41
1704 cc**********************************************************************
1705 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1706 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1707 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1708 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1710 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1711 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1712 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1713 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1715 c & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1717 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1718 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1719 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1720 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1721 c**********************************************************************
1722 apm3x(i)=g3x*d42-rimx*d43
1723 apm3y(i)=g3y*d42-rimy*d43
1724 apm3z(i)=g3z*d42-rimz*d43
1727 c*******************************************************************************
1729 c write(*,*) 'inside angvectors5'
1736 rimmx=bx(i-1)-bx(i-2)
1737 rimmy=by(i-1)-by(i-2)
1738 rimmz=bz(i-1)-bz(i-2)
1739 rim3x=bx(i-2)-bx(i-3)
1740 rim3y=by(i-2)-by(i-3)
1741 rim3z=bz(i-2)-bz(i-3)
1742 gmmx=rimmy*rimz-rimmz*rimy
1743 gmmy=rimmz*rimx-rimmx*rimz
1744 gmmz=rimmx*rimy-rimmy*rimx
1745 g3x=rim3y*rimmz-rim3z*rimmy
1746 g3y=rim3z*rimmx-rim3x*rimmz
1747 g3z=rim3x*rimmy-rim3y*rimmx
1749 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1750 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1751 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1756 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1757 cc**********************************************************************
1758 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1759 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1760 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1761 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1762 d_gmmcrimm3=d_gmmcrimm**3
1763 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1764 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1766 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1767 & .and.drim3.gt.0.0D0
1768 & .and.d_gmmcrimm3.gt.0.0D0) then
1769 d40=1.0D0/(dg30*dgmm)
1770 d41=g3gmm/(dg30*dgmm3)
1771 d42=1.0D0/(dg30*drim)
1772 d43=g3rim_/(dg30*drim3)
1773 sd40=1.0D0/(dg30*d_gmmcrimm)
1774 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1783 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1784 & -(gmmy*rimmz-gmmz*rimmy)*d41
1785 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1786 & -(gmmz*rimmx-gmmx*rimmz)*d41
1787 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1788 & -(gmmx*rimmy-gmmy*rimmx)*d41
1789 cc**********************************************************************
1790 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1791 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1792 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1793 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1795 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1796 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1797 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1798 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1800 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1801 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1802 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1803 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1804 cc**********************************************************************
1805 apm3x(i)=g3x*d42-rimx*d43
1806 apm3y(i)=g3y*d42-rimy*d43
1807 apm3z(i)=g3z*d42-rimz*d43
1817 gmx=rimy*riz-rimz*riy
1818 gmy=rimz*rix-rimx*riz
1819 gmz=rimx*riy-rimy*rix
1820 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1824 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1825 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1826 gmcrim_x=gmy*rimz-gmz*rimy
1827 gmcrim_y=gmz*rimx-gmx*rimz
1828 gmcrim_z=gmx*rimy-gmy*rimx
1829 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1830 d_gmcrim3=d_gmcrim**3
1831 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1832 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1834 if(dgm3.gt.0.0D0.and.
1835 & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1836 & .and.d_gmcrim3.gt.0.0D0) then
1837 d30=1.0D0/(dgm*dgmm)
1838 d31=gmmgm/(dgm3*dgmm)
1839 d32=gmmgm/(dgm*dgmm3)
1840 d33=1.0D0/(dgmm*dri)
1841 d34=gmmr/(dgmm3*dri)
1842 d35=gmmr/(dgmm*dri3)
1843 sd30=1.0D0/(d_gmcrim*dgmm)
1844 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1845 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1858 cc**********************************************************************
1859 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1860 & -(ciy*gmz-ciz*gmy)*d31
1861 & -(gmmy*rimmz-gmmz*rimmy)*d32
1862 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1863 & -(ciz*gmx-cix*gmz)*d31
1864 & -(gmmz*rimmx-gmmx*rimmz)*d32
1865 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1866 & -(cix*gmy-ciy*gmx)*d31
1867 & -(gmmx*rimmy-gmmy*rimmx)*d32
1868 cc**********************************************************************
1869 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1870 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1871 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1872 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1873 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1874 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1875 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1877 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1878 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1879 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1880 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1881 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1882 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1883 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1885 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1886 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1887 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1888 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1889 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1890 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1891 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1892 cc**********************************************************************
1893 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1894 & -(gmmy*rimmz-gmmz*rimmy)*d34
1896 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1897 & -(gmmz*rimmx-gmmx*rimmz)*d34
1899 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1900 & -(gmmx*rimmy-gmmy*rimmx)*d34
1904 c write(*,*) 'inside angvectors6'
1914 gx=riy*ripz-riz*ripy
1915 gy=riz*ripx-rix*ripz
1916 gz=rix*ripy-riy*ripx
1917 ggm=gmx*gx+gmy*gy+gmz*gz
1918 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1919 dg=sqrt(gx**2+gy**2+gz**2)
1925 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1927 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1928 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1929 if(dgm3.gt.0.0D0.and.
1930 & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1935 d23=1.0D0/(dgm*drip)
1936 d24=gmrp/(dgm3*drip)
1937 sd20=1.0D0/(d_gcr*dgm)
1938 sd21=gcr__gm/(d_gcr3*dgm)
1939 sd22=gcr__gm/(d_gcr*dgm3)
1950 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1951 & -(ciy*gmz-ciz*gmy)*d21
1952 & +(ripy*gz-ripz*gy)*d22
1953 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1954 & -(ciz*gmx-cix*gmz)*d21
1955 & +(ripz*gx-ripx*gz)*d22
1956 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1957 & -(cix*gmy-ciy*gmx)*d21
1958 & +(ripx*gy-ripy*gx)*d22
1959 cc**********************************************************************
1960 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1961 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1962 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1963 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1964 & +gcr_z*(-ripz*rix+gy))
1965 & -sd22*(-gmy*ciz+gmz*ciy)
1967 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1968 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1970 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1971 & -gcr_z*(ripz*riy+gx))
1972 & -sd22*(gmx*ciz-gmz*cix)
1974 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1975 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1977 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1978 & +gcr_z*(ripy*riy+ripx*rix))
1979 & -sd22*(-gmx*ciy+gmy*cix)
1980 cc**********************************************************************
1982 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1983 & -(ciy*gmz-ciz*gmy)*d24
1984 apmy(i)=(ciz*ripx-ripz*cix)*d23
1985 & -(ciz*gmx-cix*gmz)*d24
1986 apmz(i)=(cix*ripy-ripx*ciy)*d23
1987 & -(cix*gmy-ciy*gmx)*d24
1995 c-------------------------------------------------------------------------------
1996 C---------------------------------------------------------------------------------
1997 subroutine sheetforce(nca,wshet)
2000 c this should be matched with dfa.fcm
2002 parameter(maxca=800)
2003 cc**********************************************************************
2006 integer inb,nmax,iselect
2008 c real*8 dfaexp(15001)
2010 real*8 vbeta,vbetp,vbetm
2011 real*8 shefx(maxca,12)
2012 real*8 shefy(maxca,12),shefz(maxca,12)
2013 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2014 real*8 vbet(maxca,maxca)
2015 real*8 wshet(maxca,maxca)
2016 real*8 bx(maxca),by(maxca),bz(maxca)
2018 common /sheca/ bx,by,bz
2019 common /phys1/ inb,nmax,iselect
2020 common /shef/ shefx,shefy,shefz
2021 common /shee/ vbeta,vbet,vbetp,vbetm
2022 common /shetf/ shetfx,shetfy,shetfz
2039 call sheetene(nca,wshet)
2042 887 format(a,1x,i6,3x,f12.8)
2043 888 format(a,1x,i4,1x,i4,3x,f12.8)
2044 889 format(a,1x,i4,3x,f12.8)
2045 !write(2,*) 'coord : '
2047 !write(2,887) 'bx:',i,bx(i)
2048 !write(2,887) 'by:',i,by(i)
2049 !write(2,887) 'bz:',i,bz(i)
2051 !write(2,*) 'After sheetforce1'
2054 !write(2,888) 'shefx :',i,k,shefx(i,k)
2055 !write(2,888) 'shefy :',i,k,shefy(i,k)
2056 !write(2,888) 'shefz :',i,k,shefz(i,k)
2062 !write(2,*) 'After sheetforce5'
2065 !write(2,888) 'shefx :',i,k,shefx(i,k)
2066 !write(2,888) 'shefy :',i,k,shefy(i,k)
2067 !write(2,888) 'shefz :',i,k,shefz(i,k)
2073 !write(2,*) 'After sheetforce6'
2076 !write(2,888) 'shefx :',i,k,shefx(i,k)
2077 !write(2,888) 'shefy :',i,k,shefy(i,k)
2078 !write(2,888) 'shefz :',i,k,shefz(i,k)
2084 !write(2,*) 'After sheetforce11'
2087 !write(2,888) 'shefx :',i,k,shefx(i,k)
2088 !write(2,888) 'shefy :',i,k,shefy(i,k)
2089 !write(2,888) 'shefz :',i,k,shefz(i,k)
2095 !write(2,*) 'After sheetforce12'
2098 !write(2,888) 'shefx :',i,k,shefx(i,k)
2099 !write(2,888) 'shefy :',i,k,shefy(i,k)
2100 !write(2,888) 'shefz :',i,k,shefz(i,k)
2106 shetfx(i)=shetfx(i)+shefx(i,k)
2107 shetfy(i)=shetfy(i)+shefy(i,k)
2108 shetfz(i)=shetfz(i)+shefz(i,k)
2111 !write(2,*) 'Beta Finished'
2113 !write(2,889) 'shetfx : ',i,shetfx(i)
2114 !write(2,889) 'shetfy : ',i,shetfy(i)
2115 !write(2,889) 'shetfz : ',i,shetfz(i)
2121 c-------------------------------------------------------------------------------
2122 subroutine sheetene(nca,wshet)
2125 parameter(maxca=800)
2126 real*8 dfa_cutoff,dfa_cutoff_delta
2127 parameter(dfa_cutoff=15.5d0)
2128 parameter(dfa_cutoff_delta=0.5d0)
2129 cc******************************************************************************
2131 c real*8 dfaexp(15001)
2132 real*8 dtmp1, dtmp2, dtmp3
2134 real*8 vbet(maxca,maxca)
2135 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2136 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2137 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2138 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2139 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2140 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2141 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2142 real*8 cph(maxca),cth(maxca)
2143 real*8 rx(maxca,maxca)
2144 real*8 ry(maxca,maxca),rz(maxca,maxca)
2145 real*8 bx(maxca),by(maxca),bz(maxca)
2146 real*8 dis(maxca,maxca)
2148 cc**********************************************************************
2149 real*8 astx(maxca),asty(maxca),astz(maxca)
2150 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2151 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2152 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2154 real*8 wshet(maxca,maxca)
2155 real*8 dp45, dm45, w_beta
2156 real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2158 integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2160 real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2162 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2163 common /shetf/ shetfx,shetfy,shetfz
2165 common /sheca/ bx,by,bz
2166 common /phys1/ inb,nmax,iselect
2168 common /difvec/ rx,ry,rz
2169 common /coscos/ cph,cth
2170 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2171 & c00,s00,ulnex,dnex
2172 common /sheconst/ dp45,dm45,w_beta
2173 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2174 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2175 common /shee/ vbeta,vbet,vbetp,vbetm
2176 common /ulang/ ulcos
2177 cc**********************************************************************
2178 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2179 & astmmz,astm3x,astm3y,astm3z
2182 real*8 r_pair_mat(maxca,maxca)
2183 real*8 e_gcont,fprim_gcont,de_gcont
2184 ci integer istrand(maxca,maxca)
2185 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2186 ci common /shetest/ istrand,istrand_p,istrand_m
2187 common /beta_p/ r_pair_mat
2188 C-------------------------------------------------------------------------------
2192 r_pair_mat(i,j)=wshet(i,j)
2193 c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2205 if (dis(i,j).lt.dfa_cutoff) then
2206 call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2207 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2214 cc**********************************************************************
2215 y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2216 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2218 y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2219 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2232 ci if(istrand(i,j).eq.1) then
2234 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2235 yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2238 pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2239 $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2240 pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2241 $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2242 pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2243 $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2244 pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2245 $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2247 yshe1=pin1(i,j)**2+pin2(i,j)**2
2248 yshe1=-0.5d0*yshe1/dshe
2249 yshe2=pin3(i,j)**2+pin4(i,j)**2
2250 yshe2=-0.5d0*yshe2/dshe
2252 ci if((yshe1+yshe2).ge.-4) then
2259 C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2260 C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2261 C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2262 C write(*,*) 'dis(i,j):',i,j,dis(i,j)
2263 C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2264 C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2265 C write(*,*) 'pin1:',pin1(i,j)
2266 C write(*,*) 'pin2:',pin2(i,j)
2267 C write(*,*) 'pin3:',pin3(i,j)
2268 C write(*,*) 'pin4:',pin4(i,j)
2271 C write(*,*) 'yy1:',yy1
2272 C write(*,*) 'yy2:',yy2
2273 C write(*,*) 'yshe1:',yshe1
2274 C write(*,*) 'yshe2:',yshe2
2277 ci if (istrand_p(i,j).eq.1) then
2285 dtmp3 = y+yy1+yy2+yshe1+yshe2
2287 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2288 C write(*,*)'2', y,yy1,yy2
2289 C write(*,*)'3', yshe1,yshe2
2291 cc if (dtmp3.le.-35.0d0) then
2292 c vbetap(i,j)=-dp45*exp(dtmp3)
2293 cc vbetap(i,j)=0.0d0
2295 c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2296 vbetap(i,j)=-dp45*exp(dtmp3)
2299 cc if (dtmp1.le.-35.0d0) then
2300 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2301 cc vbetap1(i,j)=0.0d0
2303 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2304 c $ *dfaexp(idint(-dtmp1*1000)+1)
2305 vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2308 cc if (dtmp2.le.-35.0d0) then
2309 C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2310 cc vbetap2(i,j)=0.0d0
2312 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2313 c $ *dfaexp(idint(-dtmp2*1000)+1)
2314 vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2317 c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2318 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2319 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2321 ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2322 ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2324 ci elseif (istrand_p(i,j).eq.0)then
2330 yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2331 yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2333 pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2334 $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2335 pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2336 $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2337 pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2338 $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2339 pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2340 $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2342 yshe1=pina1(i,j)**2+pina2(i,j)**2
2343 yshe1=-0.5d0*yshe1/dshe
2344 yshe2=pina3(i,j)**2+pina4(i,j)**2
2345 yshe2=-0.5d0*yshe2/dshe
2347 ci if((yshe1+yshe2).ge.-4) then
2354 C write(*,*) 'pina1:',pina1(i,j)
2355 C write(*,*) 'pina2:',pina2(i,j)
2356 C write(*,*) 'pina3:',pina3(i,j)
2357 C write(*,*) 'pina4:',pina4(i,j)
2358 C write(*,*) 'yshe1:',yshe1
2359 C write(*,*) 'yshe2:',yshe2
2360 C write(*,*) 'dshe:',dshe
2362 ci if (istrand_m(i,j).eq.1) then
2369 dtmp3=y+yy1+yy2+yshe1+yshe2
2373 cc if(dtmp3 .le. -35.0d0) then
2374 c vbetam(i,j)=-dm45*exp(dtmp3)
2375 cc vbetam(i,j)=0.0d0
2377 c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2378 vbetam(i,j)=-dm45*exp(dtmp3)
2381 cc if(dtmp1 .le. -35.0d0) then
2382 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2383 cc vbetam1(i,j)=0.0d0
2385 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2386 c $ *dfaexp(idint(-dtmp1*1000)+1)
2387 vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2390 cc if(dtmp2.le.-35.0d0) then
2391 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2392 cc vbetam2(i,j)=0.0d0
2394 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2395 c $ *dfaexp(idint(-dtmp2*1000)+1)
2396 vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2399 ci elseif (istrand_m(i,j).eq.0)then
2406 c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2407 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2408 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2410 ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2411 ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2413 uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2414 uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2416 c write(*,*) 'uup,uum:', uup, uum
2418 c uup=vbetap1(i,j)+vbetap2(i,j)
2419 c uum=vbetam1(i,j)+vbetam2(i,j)
2424 vbeta=vbeta+vbet(i,j)*e_gcont
2427 if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
2428 c gradient correction from gcont
2429 de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
2430 shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
2431 shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
2432 shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
2434 shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
2435 shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
2436 shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
2438 c energy correction from gcont
2439 vbet(i,j)=vbet(i,j)*e_gcont
2440 vbetap(i,j)=vbetap(i,j)*e_gcont
2441 vbetap1(i,j)=vbetap1(i,j)*e_gcont
2442 vbetap2(i,j)=vbetap2(i,j)*e_gcont
2443 vbetam(i,j)=vbetam(i,j)*e_gcont
2444 vbetam1(i,j)=vbetam1(i,j)*e_gcont
2445 vbetam2(i,j)=vbetam2(i,j)*e_gcont
2449 ci elseif(istrand(i,j).eq.0)then
2453 c write(*,*) 'uup,uum:',uup,uum
2454 c write(*,*) 'vbetap(i,j):',vbetap(i,j)
2455 c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2456 c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2457 c write(*,*) 'vbetam(i,j):',vbetam(i,j)
2458 c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2459 c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2460 c write(*,*) 'uup:',uup
2461 c write(*,*) 'uum:',uum
2462 c write(*,*) 'vbetp:',vbetp
2463 c write(*,*) 'vbetm:',vbetm
2464 c write(*,*) 'vbet(i,j):',vbet(i,j)
2481 ! write(*,*) 'I,J:', i,j
2482 ! write(*,*) 'vbetap(i,j):',vbetap(i,j)
2483 ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2484 ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2485 ! write(*,*) 'vbetam(i,j):',vbetam(i,j)
2486 ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2487 ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2488 ! write(*,*) 'vbet(i,j):',vbet(i,j)
2494 c-------------------------------------------------------------------------------
2495 subroutine sheetforce1
2498 parameter(maxca=800)
2499 real*8 dfa_cutoff,dfa_cutoff_delta
2500 parameter(dfa_cutoff=15.5d0)
2501 parameter(dfa_cutoff_delta=0.5d0)
2502 cc**********************************************************************
2503 real*8 vbet(maxca,maxca)
2504 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2505 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2506 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2507 real*8 cph(maxca),cth(maxca)
2508 real*8 rx(maxca,maxca)
2509 real*8 ry(maxca,maxca),rz(maxca,maxca)
2510 real*8 bx(maxca),by(maxca),bz(maxca)
2511 real*8 dis(maxca,maxca)
2512 real*8 shefx(maxca,12)
2513 real*8 shefy(maxca,12),shefz(maxca,12)
2514 real*8 atx(maxca),aty(maxca),atz(maxca)
2515 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2516 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2517 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2518 real*8 apx(maxca),apy(maxca),apz(maxca)
2519 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2520 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2521 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2523 real*8 astx(maxca),asty(maxca),astz(maxca)
2524 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2525 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2526 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2528 real*8 w_beta,dp45, dm45
2529 real*8 vbeta, vbetp, vbetm
2530 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2531 $ c00,s00,ulnex,dnex
2532 integer inb,nmax,iselect
2534 common /phys1/ inb,nmax,iselect
2536 common /difvec/ rx,ry,rz
2537 common /coscos/ cph,cth
2538 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2539 $ c00,s00,ulnex,dnex
2540 common /sheconst/ dp45,dm45,w_beta
2541 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2542 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2543 $ atmmz,atm3x,atm3y,atm3z
2544 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2545 $ apmmz,apm3x,apm3y,apm3z
2546 common /shef/ shefx,shefy,shefz
2547 common /shee/ vbeta,vbet,vbetp,vbetm
2548 common /ulang/ ulcos
2549 c c**********************************************************************
2550 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2551 $ astmmz,astm3x,astm3y,astm3z
2553 C--------------------------------------------------------------------------------
2555 integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2556 real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2557 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2558 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2559 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2560 real*8 e_gcont,fprim_gcont
2561 C--------------------------------------------------------------------------------
2566 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2571 cc1=(ulcos(imm)-ulnex)/dnex
2572 dmm=cc1/(dis(imm,im)*dis(im,i))
2573 dmm__=cc1*ulcos(imm)/dis(im,i)**2
2574 fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2575 fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2576 fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2580 fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2581 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2582 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2592 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2597 cc1=(ulcos(imm)-ulnex)/dnex
2598 cc2=(ulcos(im)-ulnex)/dnex
2599 dmm1=cc1/(dis(imm,im)*dis(im,i))
2600 dmm2=cc2/(dis(im,i)*dis(i,ip))
2601 dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2602 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2603 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2604 cc**********************************************************************
2605 fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2606 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2607 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2608 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2609 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2610 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2614 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2615 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2616 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2625 c3=(cth(im)*c00+sth(im)*s00-1)/dca
2630 cc2=(ulcos(im)-ulnex)/dnex
2631 cc3=(ulcos(i)-ulnex)/dnex
2632 dmm2=cc2/(dis(im,i)*dis(i,ip))
2633 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2634 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2635 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2636 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2637 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2638 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2639 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2640 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2641 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2642 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2646 fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2647 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2648 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2656 c4=(cth(i)*c00+sth(i)*s00-1)/dca
2661 cc3=(ulcos(i)-ulnex)/dnex
2662 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2663 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2664 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2665 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2666 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2670 fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2671 fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2672 fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2681 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2686 cc7=(ulcos(jmm)-ulnex)/dnex
2687 dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2688 dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2689 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2690 fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2691 fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2695 fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2696 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2697 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2706 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2711 cc7=(ulcos(jmm)-ulnex)/dnex
2712 cc8=(ulcos(jm)-ulnex)/dnex
2713 dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2714 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2715 dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2716 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2717 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2718 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2719 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2720 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2721 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2722 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2723 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2727 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2728 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2729 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2739 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2744 cc8=(ulcos(jm)-ulnex)/dnex
2745 cc9=(ulcos(j)-ulnex)/dnex
2746 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2747 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2748 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2749 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2750 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2751 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2752 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2753 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2754 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2755 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2756 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2760 fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2761 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2762 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2771 c10=(cth(j)*c00+sth(j)*s00-1)/dca
2776 cc9=(ulcos(j)-ulnex)/dnex
2777 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2778 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2779 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2780 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2781 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2785 fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2786 fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2787 fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2795 c----------------------------------------------------------------------------
2796 subroutine sheetforce5
2799 parameter(maxca=800)
2800 real*8 dfa_cutoff,dfa_cutoff_delta
2801 parameter(dfa_cutoff=15.5d0)
2802 parameter(dfa_cutoff_delta=0.5d0)
2803 cc**********************************************************************
2804 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2805 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2806 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2807 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2808 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2809 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2810 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2811 real*8 rx(maxca,maxca)
2812 real*8 ry(maxca,maxca),rz(maxca,maxca)
2813 real*8 bx(maxca),by(maxca),bz(maxca)
2814 real*8 dis(maxca,maxca)
2815 real*8 shefx(maxca,12),shefy(maxca,12)
2816 real*8 shefz(maxca,12)
2817 real*8 dp45,dm45,w_beta
2818 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2819 $ c00,s00,ulnex,dnex
2820 integer inb,nmax,iselect
2821 cc**********************************************************************
2822 common /phys1/ inb,nmax,iselect
2824 common /difvec/ rx,ry,rz
2825 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2826 $ c00,s00,ulnex,dnex
2827 common /sheconst/ dp45,dm45,w_beta
2828 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2829 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2830 common /shef/ shefx,shefy,shefz
2831 ci integer istrand(maxca,maxca)
2832 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2833 ci common /shetest/ istrand,istrand_p,istrand_m
2834 c********************************************************************************
2836 integer i,imm,im,jp,jpp,j
2837 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2838 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2839 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2840 real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2841 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2842 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2843 real*8 e_gcont,fprim_gcont
2844 c********************************************************************************
2850 if (dis(imm,j).lt.dfa_cutoff) then
2851 call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2852 & dfa_cutoff_delta,e_gcont,fprim_gcont)
2857 ci if(istrand(imm,j).eq.1
2858 ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2861 yy1=-(dis(i,jpp)-ulhb)/dlhb
2862 y1x=rx(jpp,i)/dis(i,jpp)
2863 y1y=ry(jpp,i)/dis(i,jpp)
2864 y1z=rz(jpp,i)/dis(i,jpp)
2869 yy33=1.0D0/(dis(im,jp)*dis(im,i))
2870 yyy3=pin1(imm,j)/(dis(im,i)**2)
2871 yy3=-pin1(imm,j)/dshe
2872 y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2873 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2874 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2876 yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2877 yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2878 yyy4b=pin3(imm,j)/(dis(im,i)**2)
2879 yy4=-pin3(imm,j)/dshe
2880 y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2881 $ -yyy4b*rx(im,i))*yy4
2882 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2883 $ -yyy4b*ry(im,i))*yy4
2884 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2885 $ -yyy4b*rz(im,i))*yy4
2888 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2889 yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2890 yy5=-pin4(imm,j)/dshe
2891 y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2892 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2893 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2906 shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2907 $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2908 shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2909 $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2910 shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2911 $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2913 ! shefx(i,5)=shefx(i,5)
2914 ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2915 ! shefy(i,5)=shefy(i,5)
2916 ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2917 ! shefz(i,5)=shefz(i,5)
2918 ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2920 yy6=-(dis(i,jp)-uldhb)/dldhb
2921 y6x=rx(jp,i)/dis(i,jp)
2922 y6y=ry(jp,i)/dis(i,jp)
2923 y6z=rz(jp,i)/dis(i,jp)
2928 yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2929 yyy8=pina1(imm,j)/(dis(im,i)**2)
2930 yy8=-pina1(imm,j)/dshe
2931 y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2932 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2933 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2935 yy99=1.0D0/(dis(jp,i)*dis(im,i))
2936 yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2937 yyy9b=pina3(imm,j)/(dis(im,i)**2)
2938 yy9=-pina3(imm,j)/dshe
2939 y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2940 $ -yyy9b*rx(im,i))*yy9
2941 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2942 $ -yyy9b*ry(im,i))*yy9
2943 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2944 $ -yyy9b*rz(im,i))*yy9
2946 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2947 yyy10=pina4(imm,j)/(dis(jp,i)**2)
2948 yy10=-pina4(imm,j)/dshe
2949 y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2950 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2951 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2953 sx=y66x+y8x+y9x+y10x
2954 sy=y66y+y8y+y9y+y10y
2955 sz=y66z+y8z+y9z+y10z
2964 shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
2965 $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2966 shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
2967 $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2968 shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
2969 $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2971 ! shefx(i,5)=shefx(i,5)
2972 ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2973 ! shefy(i,5)=shefy(i,5)
2974 ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2975 ! shefz(i,5)=shefz(i,5)
2976 ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2986 c--------------------------------------------------------------------------c
2987 subroutine sheetforce6
2990 parameter(maxca=800)
2991 real*8 dfa_cutoff,dfa_cutoff_delta
2992 parameter(dfa_cutoff=15.5d0)
2993 parameter(dfa_cutoff_delta=0.5d0)
2994 cc**********************************************************************
2995 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2996 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2997 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2998 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2999 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3000 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3001 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3002 real*8 rx(maxca,maxca)
3003 real*8 ry(maxca,maxca),rz(maxca,maxca)
3004 real*8 bx(maxca),by(maxca),bz(maxca)
3005 real*8 dis(maxca,maxca)
3006 real*8 shefx(maxca,12),shefy(maxca,12)
3007 real*8 shefz(maxca,12)
3008 real*8 dp45,dm45,w_beta
3009 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3010 $ c00,s00,ulnex,dnex
3011 integer inb,nmax,iselect
3012 cc**********************************************************************
3013 common /phys1/ inb,nmax,iselect
3015 common /difvec/ rx,ry,rz
3016 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3017 $ c00,s00,ulnex,dnex
3018 common /sheconst/ dp45,dm45,w_beta
3019 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3020 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3021 common /shef/ shefx,shefy,shefz
3022 ci integer istrand(maxca,maxca)
3023 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3024 ci common /shetest/ istrand,istrand_p,istrand_m
3025 cc**********************************************************************
3027 integer i,imm,im,jp,jpp,j,ip
3028 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3029 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3030 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3031 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3032 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
3033 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
3034 real*8 e_gcont,fprim_gcont
3035 C********************************************************************************
3041 if (dis(im,j).lt.dfa_cutoff) then
3042 call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3043 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3048 ci if(istrand(im,j).eq.1
3049 ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
3052 yy1=-(dis(i,jp)-ulhb)/dlhb
3053 y1x=rx(jp,i)/dis(i,jp)
3054 y1y=ry(jp,i)/dis(i,jp)
3055 y1z=rz(jp,i)/dis(i,jp)
3060 yy33=1.0D0/(dis(i,jp)*dis(i,ip))
3061 yyy3a=pin1(im,j)/(dis(i,jp)**2)
3062 yyy3b=pin1(im,j)/(dis(i,ip)**2)
3063 yy3=-pin1(im,j)/dshe
3064 y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
3065 $ +yyy3b*rx(i,ip))*yy3
3066 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
3067 $ +yyy3b*ry(i,ip))*yy3
3068 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
3069 $ +yyy3b*rz(i,ip))*yy3
3071 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3072 yyy4=pin2(im,j)/(dis(i,jp)**2)
3073 yy4=-pin2(im,j)/dshe
3074 y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3075 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3076 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3078 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3079 yyy5=pin3(im,j)/(dis(i,ip)**2)
3080 yy5=-pin3(im,j)/dshe
3081 y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3082 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3083 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3096 shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3097 $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3098 shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3099 $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3100 shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3101 $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3102 ! shefx(i,6)=shefx(i,6)
3103 ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3104 ! shefy(i,6)=shefy(i,6)
3105 ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3106 ! shefz(i,6)=shefz(i,6)
3107 ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3109 yy6=-(dis(jpp,i)-uldhb)/dldhb
3110 y6x=rx(jpp,i)/dis(jpp,i)
3111 y6y=ry(jpp,i)/dis(jpp,i)
3112 y6z=rz(jpp,i)/dis(jpp,i)
3117 yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3118 yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3119 yyy8b=pina1(im,j)/(dis(i,ip)**2)
3120 yy8=-pina1(im,j)/dshe
3121 y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3122 $ +yyy8b*rx(i,ip))*yy8
3123 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3124 $ +yyy8b*ry(i,ip))*yy8
3125 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3126 $ +yyy8b*rz(i,ip))*yy8
3128 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3129 yyy9=pina2(im,j)/(dis(i,jpp)**2)
3130 yy9=-pina2(im,j)/dshe
3131 y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3132 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3133 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3135 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3136 yyy10=pina3(im,j)/(dis(i,ip)**2)
3137 yy10=-pina3(im,j)/dshe
3138 y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3139 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3140 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3142 sx=y66x+y8x+y9x+y10x
3143 sy=y66y+y8y+y9y+y10y
3144 sz=y66z+y8z+y9z+y10z
3153 shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3154 $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3155 shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3156 $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3157 shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3158 $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3160 ! shefx(i,6)=shefx(i,6)
3161 ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3162 ! shefy(i,6)=shefy(i,6)
3163 ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3164 ! shefz(i,6)=shefz(i,6)
3165 ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3175 c-----------------------------------------------------------------------
3176 subroutine sheetforce11
3179 parameter(maxca=800)
3180 real*8 dfa_cutoff,dfa_cutoff_delta
3181 parameter(dfa_cutoff=15.5d0)
3182 parameter(dfa_cutoff_delta=0.5d0)
3183 cc**********************************************************************
3184 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3185 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3186 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3187 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3188 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3189 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3190 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3191 real*8 rx(maxca,maxca)
3192 real*8 ry(maxca,maxca),rz(maxca,maxca)
3193 real*8 bx(maxca),by(maxca),bz(maxca)
3194 real*8 dis(maxca,maxca)
3195 real*8 shefx(maxca,12),shefy(maxca,12)
3196 real*8 shefz(maxca,12)
3197 real*8 dp45,dm45,w_beta
3198 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3199 $ c00,s00,ulnex,dnex
3200 integer inb,nmax,iselect
3201 cc**********************************************************************
3202 common /phys1/ inb,nmax,iselect
3204 common /difvec/ rx,ry,rz
3205 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3206 $ c00,s00,ulnex,dnex
3207 common /sheconst/ dp45,dm45,w_beta
3208 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3209 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3210 common /shef/ shefx,shefy,shefz
3211 ci integer istrand(maxca,maxca)
3212 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3213 ci common /shetest/ istrand,istrand_p,istrand_m
3214 C********************************************************************************
3216 integer j,jm,jmm,ip,i,ipp
3217 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3218 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3219 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3220 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3221 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3222 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3223 real*8 e_gcont,fprim_gcont
3224 C********************************************************************************
3231 if (dis(i,jmm).lt.dfa_cutoff) then
3232 call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3233 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3238 ci if(istrand(i,jmm).eq.1
3239 ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3242 yy1=-(dis(ipp,j)-ulhb)/dlhb
3243 y1x=rx(ipp,j)/dis(ipp,j)
3244 y1y=ry(ipp,j)/dis(ipp,j)
3245 y1z=rz(ipp,j)/dis(ipp,j)
3250 yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3251 yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3252 yy3=-pin2(i,jmm)/dshe
3253 y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3254 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3255 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3257 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3258 yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3259 yy4=-pin3(i,jmm)/dshe
3260 y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3261 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3262 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3264 yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3265 yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3266 yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3267 yy5=-pin4(i,jmm)/dshe
3268 y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3269 $ -yyy5b*rx(jm,j))*yy5
3270 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3271 $ -yyy5b*ry(jm,j))*yy5
3272 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3273 $ -yyy5b*rz(jm,j))*yy5
3286 shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3287 $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3288 shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3289 $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3290 shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3291 $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3293 ! shefx(j,11)=shefx(j,11)
3294 ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3295 ! shefy(j,11)=shefy(j,11)
3296 ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3297 ! shefz(j,11)=shefz(j,11)
3298 ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3300 yy6=-(dis(ip,j)-uldhb)/dldhb
3301 y6x=rx(ip,j)/dis(ip,j)
3302 y6y=ry(ip,j)/dis(ip,j)
3303 y6z=rz(ip,j)/dis(ip,j)
3308 yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3309 yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3310 yy8=-pina1(i,jmm)/dshe
3311 y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3312 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3313 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3315 yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3316 yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3317 yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3318 yy9=-pina2(i,jmm)/dshe
3319 y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3320 $ -yyy9b*rx(jm,j))*yy9
3321 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3322 $ -yyy9b*ry(jm,j))*yy9
3323 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3324 $ -yyy9b*rz(jm,j))*yy9
3326 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3327 yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3328 yy10=-pina4(i,jmm)/dshe
3329 y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3330 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3331 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3333 sx=y66x+y8x+y9x+y10x
3334 sy=y66y+y8y+y9y+y10y
3335 sz=y66z+y8z+y9z+y10z
3344 shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3345 $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3346 shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3347 $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3348 shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3349 $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3351 ! shefx(j,11)=shefx(j,11)
3352 ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3353 ! shefy(j,11)=shefy(j,11)
3354 ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3355 ! shefz(j,11)=shefz(j,11)
3356 ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3366 c-----------------------------------------------------------------------
3367 subroutine sheetforce12
3370 parameter(maxca=800)
3371 real*8 dfa_cutoff,dfa_cutoff_delta
3372 parameter(dfa_cutoff=15.5d0)
3373 parameter(dfa_cutoff_delta=0.5d0)
3374 cc**********************************************************************
3375 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3376 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3377 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3378 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3379 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3380 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3381 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3382 real*8 rx(maxca,maxca)
3383 real*8 ry(maxca,maxca),rz(maxca,maxca)
3384 real*8 bx(maxca),by(maxca),bz(maxca)
3385 real*8 dis(maxca,maxca)
3386 real*8 shefx(maxca,12),shefy(maxca,12)
3387 real*8 shefz(maxca,12)
3388 real*8 dp45,dm45,w_beta
3389 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3390 $ c00,s00,ulnex,dnex
3391 integer inb,nmax,iselect
3392 cc**********************************************************************
3393 common /phys1/ inb,nmax,iselect
3395 common /difvec/ rx,ry,rz
3396 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3397 $ c00,s00,ulnex,dnex
3398 common /sheconst/ dp45,dm45,w_beta
3399 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3400 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3401 common /shef/ shefx,shefy,shefz
3402 ci integer istrand(maxca,maxca)
3403 ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3404 ci common /shetest/ istrand,istrand_p,istrand_m
3405 cc**********************************************************************
3407 integer j,jm,jmm,ip,i,ipp,jp
3408 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3409 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3410 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3411 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3412 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3413 real*8 e_gcont,fprim_gcont
3414 !c*************************************************************************c
3420 if (dis(i,jm).lt.dfa_cutoff) then
3421 call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3422 & dfa_cutoff_delta,e_gcont,fprim_gcont)
3427 ci if(istrand(i,jm).eq.1
3428 ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3431 yy1=-(dis(ip,j)-ulhb)/dlhb
3432 y1x=rx(ip,j)/dis(ip,j)
3433 y1y=ry(ip,j)/dis(ip,j)
3434 y1z=rz(ip,j)/dis(ip,j)
3439 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3440 yyy3=pin1(i,jm)/(dis(ip,j)**2)
3441 yy3=-pin1(i,jm)/dshe
3442 y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3443 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3444 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3445 yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3447 yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3448 yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3449 yy4=-pin2(i,jm)/dshe
3450 y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3451 $ +yyy4b*rx(j,jp))*yy4
3452 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3453 $ +yyy4b*ry(j,jp))*yy4
3454 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3455 $ +yyy4b*rz(j,jp))*yy4
3457 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3458 yyy5=pin4(i,jm)/(dis(j,jp)**2)
3459 yy5=-pin4(i,jm)/dshe
3460 y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3461 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3462 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3475 shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3476 $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3477 shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3478 $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3479 shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3480 $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3482 ! shefx(j,12)=shefx(j,12)
3483 ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3484 ! shefy(j,12)=shefy(j,12)
3485 ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3486 ! shefz(j,12)=shefz(j,12)
3487 ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3489 yy6=-(dis(ipp,j)-uldhb)/dldhb
3490 y6x=rx(ipp,j)/dis(ipp,j)
3491 y6y=ry(ipp,j)/dis(ipp,j)
3492 y6z=rz(ipp,j)/dis(ipp,j)
3497 yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3498 yyy8=pina2(i,jm)/(dis(j,jp)**2)
3499 yy8=-pina2(i,jm)/dshe
3500 y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3501 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3502 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3504 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3505 yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3506 yy9=-pina3(i,jm)/dshe
3507 y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3508 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3509 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3511 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3512 yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3513 yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3514 yy10=-pina4(i,jm)/dshe
3515 y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3516 $ +yyy10b*rx(j,jp))*yy10
3517 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3518 $ +yyy10b*ry(j,jp))*yy10
3519 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3520 $ +yyy10b*rz(j,jp))*yy10
3522 sx=y66x+y8x+y9x+y10x
3523 sy=y66y+y8y+y9y+y10y
3524 sz=y66z+y8z+y9z+y10z
3533 shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3534 $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3535 shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3536 $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3537 shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3538 $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3549 C===============================================================================