1 double precision function rdif(n,x,g,ider)
4 include "DIMENSIONS.ZSCOPT"
6 parameter (nmax=max_paropt)
11 double precision g_(max_paropt),sum_
13 include "COMMON.IOUNITS"
14 include "COMMON.PMFDERIV"
15 include "COMMON.WEIGHTS"
16 include "COMMON.TORSION"
17 double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
18 & b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0:3),
20 & eello_turn3_b2(3,2),eello_turn3_b1(3,2),eello_turn3_e1(2,2,2),
21 & eello_turn3_e2(2,2,2),eello_turn3_e0m1(3),eello_turn3_e0m2(3),
22 & c2(0:3),em(2,2,2),em2(2,2,2),e0m(3),e0m2(3),jac(nmax)
23 double precision v1_kcc_loc(3,3,2),v2_kcc_loc(3,3,2),
24 & etoriv1(3,3,2),etoriv2(3,3,2),d1,d2
25 double precision eneelloc3(4),eneelloc3b(3,2,2,2,4)
26 double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
27 & phi,cost1,cost2,sint1,sint2,
28 & cosg,sing,cos2g,sin2g,sum,d3,delta,f,fb11,fb21,fb12,fb22,fc11,
29 & fc12,fd11,fd12,etori,sint1sint2,cosphi,sinphi,sumvalc,sumvals
30 double precision conv /0.01745329251994329576d0/
31 double precision r0pp(3) /5.2739d0,5.4561d0,5.2261d0/
32 double precision epspp(3) /1.6027d0,1.4879d0,0.0779d0/
33 double precision dCACA(2) /3.784904d0,3.810180d0/
34 double precision facp1(4),facpp(4),facturn3,ddist
36 double precision x(n),g(n)
37 integer ipoint,iipoint,ip,it1,it2,itp1,itp2,i,j,k,l,ii,jj!,irt1,irt2
40 c nfun.eq.-1 : no broadcast
41 c n.eq.-1 : stop function evaluation
42 c nfun.ne.-1 .and. n.ne.-1 : master distribute variables to slaves, all do their share
44 c if (nfun.ne.-1) then
45 c call MPI_Bcast(n,1,MPI_INTEGER,Master,Comm1,ierror)
48 call MPI_Bcast(ider,1,MPI_INTEGER,Master,Comm1,ierror)
50 c call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION,Master,
53 c write (2,*) "x",(x(i),i=1,n)
56 c write (iout,*) "ider",ider," lder",lder
58 c print *,"rdif npoint",npoint
60 c Zero out gradient and hessian
73 c-- Following for ala-ala only
86 do it2=-nloctyp+1,nloctyp-1
88 write (iout,*) "it1",it1," it2",it2," e0",e0(it1,it2)
90 c following to test with selected types
94 if (it1.lt.nloctyp-1) then
99 if (iabs(it2).lt.nloctyp-1) then
108 facturn3=wturn_PMF*dsqrt(epspp(1))*4.2d0**3
110 facturn3=wturn_PMF*dsqrt(epspp(2))*4.2d0**3
118 facp1(4)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
120 facp1(2)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
122 facp1(2)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
125 facp1(3)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
127 facp1(3)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
129 if (itp1+itp2.eq.2) then
130 facp1(1)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
131 else if (itp1+itp2.eq.3) then
132 facp1(1)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
134 facp1(1)=wello_PMF*dsqrt(epspp(3))*4.2d0**3
136 c print *,"itp1",itp1," itp2",itp2," facp1",facp1
139 c irt2=15+20*iabs(it2)
140 c-- The following for ALL residues
142 c irt2=17+31*iabs(it2)
143 c-- The followinig for ala-ala only
146 c print *,"it1",it1," it2",it2
147 c print *,irt1,irt2,irt1+31,irt2+31
148 c----------------------
149 c if (FOURIERSPLIT) then
152 b1m(k,i)=bnew1tor(k,i,it2)
153 b1m2(k,i)=bnew2tor(k,i,it2)
154 b2m1(k,i)=bnew1tor(k,i,it1)
155 b2m(k,i)=bnew2tor(k,i,it1)
160 cm(k,i)=ccnewtor(k,i,it2)
161 dm(k,i)=ddnewtor(k,i,it1)
167 em(k,j,i)=eenewtor(k,j,i,it1)
168 em2(k,j,i)=eenewtor(k,j,i,it2)
173 e0m(k)=e0newtor(k,it1)
174 e0m2(k)=e0newtor(k,it2)
176 c Invert coefficients for L-chirality
178 write(iout,*) "it1",it1," it2",it2
179 write(iout,*) "b1i (b1m)"
181 write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
183 write(iout,*) "b2j (b1m2)"
185 write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
187 write(iout,*) "b1j (b2m1)"
189 write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
191 write(iout,*) "b2j (b2m)"
193 write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
197 write (iout,'(2f10.5)') (cm(i,j),j=1,2)
201 write (iout,'(2f10.5)') (dm(i,j),j=1,2)
203 print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
213 v1_kcc_loc(k,l,1)=b2m(l,1)*b1m(k,1)+b2m(l,2)*b1m(k,2)
214 v2_kcc_loc(k,l,1)=b2m(l,2)*b1m(k,1)+b2m(l,1)*b1m(k,2)
219 v1_kcc_loc(k,l,2)=-(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
220 v2_kcc_loc(k,l,2)=-(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
225 if (torsion_pmf) then
227 c write (iout,*) "it1",it1," it2",it2," etori",e0(it1,it2),
228 c & " mask",mask_e0(it1,it2)
229 c write (iout,*)"mask_bnew1 1",it2,
230 c & (mask_bnew1tor(j,1,iabs(it2)),j=1,3)
231 c write (iout,*)"mask_bnew1 2",it2,
232 c & (mask_bnew1tor(j,2,iabs(it2)),j=1,3)
233 c write (iout,*)"mask_bnew2 1",it1,(mask_bnew2tor(j,1,it1),j=1,3)
234 c write (iout,*)"mask_bnew2 2",it1,(mask_bnew1tor(j,2,it1),j=1,3)
241 gam=(zz(3,i)+delta-180)*conv
251 b11=(b1m(1,1)+b1m(2,1)*cost2+b1m(3,1)*cost2*cost2)*sint2
252 b12=(b1m(1,2)+b1m(2,2)*cost2+b1m(3,2)*cost2*cost2)*sint2
253 b21=(b2m(1,1)+b2m(2,1)*cost1+b2m(3,1)*cost1*cost1)*sint1
254 b22=(b2m(1,2)+b2m(2,2)*cost1+b2m(3,2)*cost1*cost1)*sint1
255 d11=(dm(1,1)+dm(2,1)*cost1)*sint1*sint1
256 d12=(dm(1,2)+dm(2,2)*cost1)*sint1*sint1
257 c11=(cm(1,1)+cm(2,1)*cost2)*sint2*sint2
258 c12=(cm(1,2)+cm(2,2)*cost2)*sint2*sint2
260 write (iout,*) "b11",b11," b12",b12," b21",b21," b22",b22
261 write (iout,*) "d11",d11," d12",d12," c11",c11," c12",c12
262 write (iout,*) "cosg",cosg," cos2g",cos2g," sing",sing,
265 f=e0(it1,it2)+(b21*b11+b22*b12)*cosg+(b22*b11+b21*b12)*sing
266 & -(c11*d11-c12*d12)*cos2g-(c12*d11+c11*d12)*sin2g
267 fb11=b21*cosg+b22*sing
268 fb12=b22*cosg+b21*sing
269 fb21=b11*cosg+b12*sing
270 fb22=b12*cosg+b11*sing
271 fc11=-d11*cos2g-d12*sin2g
272 fc12= d12*cos2g-d11*sin2g
273 fd11=-c11*cos2g-c12*sin2g
274 fd12= c12*cos2g-c11*sin2g
276 c check with the formula from unres
290 sint1sint2=sint1sint2*sint1*sint2
295 sumvalc=sumvalc+v1_kcc_loc(l,k,j)*c1(k)*c2(l)
296 sumvals=sumvals+v2_kcc_loc(l,k,j)*c1(k)*c2(l)
297 etoriv1(l,k,j)=c1(k)*c2(l)*sint1sint2*cosphi
298 etoriv2(l,k,j)=c1(k)*c2(l)*sint1sint2*sinphi
301 c if (j.eq.1) print *,"j",j," sumvlac",sumvalc*sint1sint2,
303 c & "sumvals",sumvals*sint1sint2,b22*b11+b21*b12
304 c if (j.eq.2) print *,"j",j," sumvlac",sumvalc*sint1sint2,
305 c & -(c11*d11-c12*d12),
306 c & "sumvals",sumvals*sint1sint2,-(c12*d11+c11*d12)
307 etori=etori+(sumvalc*cosphi+sumvals*sinphi)*sint1sint2
311 c v1_kcc(k,l,1)=b1m(k,1)*b2m(l,1)+b1m(k,2)*b2m(l,2)
312 c v2_kcc(k,l,1)=b1m(k,1)*b2m(l,2)+b1m(k,2)*b2m(l,1)
317 c v1_kcc(k,l,2)=-0.25d0*(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
318 c v2_kcc(k,l,2)=-0.25d0*(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
321 c print *,"i",i," theta",zz(i,1),zz(i,2)," gamma",zz(i,3),
322 c & " f",f," etori",etori
327 write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
328 & y(1,i),f,yc(1,i),y(1,i)-yc(1,i)
330 diff(1,i)=y(1,i)-yc(1,i)
331 sum=sum+diff(1,i)*diff(1,i)*w(1,i)
333 c Compute derivatives of the energy in torsional parameters
341 if (mask_e0(it1,it2).gt.0) jac(mask_e0(it1,it2))=1.0d0
342 c Derivatives in b1 of the second residue
346 jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
349 jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
350 & +etoriv2(j,k,1)*b2m(k,2)
353 jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
357 jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,2)
358 & +etoriv2(j,k,1)*b2m(k,1)
359 else if (it2.lt.0) then
360 jac(jj)=jac(jj)-etoriv1(j,k,1)*b2m(k,2)
361 & -etoriv2(j,k,1)*b2m(k,1)
366 c Derivatives in b2 of the first residue
369 jj=iabs(mask_bnew2tor(j,1,it1))
372 jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
373 & +etoriv2(k,j,1)*b1m(k,2)
376 jj=iabs(mask_bnew2tor(j,2,it1))
380 jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
381 & +etoriv2(k,j,1)*b1m(k,1)
389 jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
392 c jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,1)
393 c & -0.25d0*etoriv2(j,k,2)*dm(k,2)
394 jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,1)
395 & -etoriv2(j,k,2)*dm(k,2)
398 jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
402 c jac(jj+1)=jac(jj)+0.25d0*etoriv1(j,k,2)*dm(k,2)
403 c & -0.25d0*etoriv2(j,k,2)*dm(k,1)
404 jac(jj)=jac(jj)+etoriv1(j,k,2)*dm(k,2)
405 & -etoriv2(j,k,2)*dm(k,1)
406 else if (it2.lt.0) then
407 c jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,2)
408 c & +0.25d0*etoriv2(j,k,2)*dm(k,1)
409 jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,2)
410 & +etoriv2(j,k,2)*dm(k,1)
418 jj=iabs(mask_ddnewtor(j,1,it1))
421 c jac(jj)=jac(jj)-0.25d0*etoriv1(k,j,2)*cm(k,1)
422 c & -0.25d0*etoriv2(k,j,2)*cm(k,2)
423 jac(jj)=jac(jj)-etoriv1(k,j,2)*cm(k,1)
424 & -etoriv2(k,j,2)*cm(k,2)
427 jj=iabs(mask_ddnewtor(j,2,it1))
431 c jac(jj)=jac(jj)+0.25d0*etoriv1(k,j,2)*cm(k,2)
432 c & -0.25d0*etoriv2(k,j,2)*cm(k,1)
433 jac(jj)=jac(jj)+etoriv1(k,j,2)*cm(k,2)
434 & -etoriv2(k,j,2)*cm(k,1)
439 c Contributions to the gradient and the hessian
441 g(k)=g(k)+diff(1,i)*jac(k)*w(1,i)
442 c h(k,k)=h(k,k)+jac(k)**2*w(1,i)
444 c h(k,l)=h(k,l)+jac(k)*jac(l)*w(1,i)
448 print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
459 b1m(k,i)=bnew1(k,i,it2)
460 b1m2(k,i)=bnew2(k,i,it2)
461 b2m1(k,i)=bnew1(k,i,it1)
462 b2m(k,i)=bnew2(k,i,it1)
468 em(k,j,i)=eenew(k,j,i,it1)
469 em2(k,j,i)=eenew(k,j,i,it2)
486 gam=(zz(3,i)+delta-180)*conv
489 c eturn3 contributions
490 call eturn3PMF(theta1,theta2,gam,b2m1,b1m2,em,em2,e0m,e0m2,
491 & facturn3,d1,d2,d3,eello_turn3,eello_turn3_b2,eello_turn3_b1,
492 & eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
495 write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
496 & y(2,i),yc(2,i),y(2,i)-yc(2,i)
498 diff(2,i)=y(2,i)-yc(2,i)
499 sum=sum+diff(2,i)*diff(2,i)*w(2,i)
504 if (mask_turn_PMF.gt.0) jac(mask_turn_PMF)=eello_turn3/wturn_PMF
505 c Derivatives of in b1 of the first residue
508 jj=mask_bnew1(j,1,it1)
511 jac(jj)=jac(jj)+eello_turn3_b2(j,1)
515 jj=mask_bnew1(j,2,it1)
518 jac(jj)=jac(jj)+eello_turn3_b2(j,2)
522 c Derivatives of in b2 of the second residue
524 jj=mask_bnew2(j,1,iabs(it2))
526 jac(jj)=jac(jj)+eello_turn3_b1(j,1)
530 jj=mask_bnew2(j,2,iabs(it2))
533 jac(jj)=jac(jj)+eello_turn3_b1(j,2)
534 else if (it2.lt.0) then
535 jac(jj)=jac(jj)-eello_turn3_b1(j,2)
541 c Derivatives in e of the first residue
544 jj=mask_eenew(j,1,1,it1)
545 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,1)
546 jj=mask_eenew(j,2,2,it1)
547 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,2)
549 jj=mask_eenew(j,1,2,it1)
550 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,2)
551 jj=mask_eenew(j,2,1,it1)
552 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,1)
557 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(1)
560 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(2)
562 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(3)
564 c Derivatives in e of the second residue
567 jj=mask_eenew(j,1,1,iabs(it2))
568 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,1)
569 jj=mask_eenew(j,2,2,iabs(it2))
570 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,2)
572 jj=mask_eenew(j,1,2,iabs(it2))
573 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,2)
574 jj=mask_eenew(j,2,1,iabs(it2))
575 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,1)
576 else if (it2.lt.0) then
577 jj=mask_eenew(j,1,2,iabs(it2))
578 if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,1,2)
579 jj=mask_eenew(j,2,1,iabs(it2))
580 if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,2,1)
584 jj=mask_e0new(1,iabs(it2))
585 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(1)
587 jj=mask_e0new(2,iabs(it2))
588 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(2)
589 jj=mask_e0new(3,iabs(it2))
590 if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(3)
591 else if (it2.lt.0) then
592 jj=mask_e0new(2,iabs(it2))
593 if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(2)
594 jj=mask_e0new(3,iabs(it2))
595 if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(3)
598 c Contributions to the gradient and the hessian
600 g(k)=g(k)+diff(2,i)*jac(k)*w(2,i)
601 c h(k,k)=h(k,k)+jac(k)**2*w(2,i)
603 c h(k,l)=h(k,l)+jac(k)*jac(l)*w(2,i)
607 print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
622 do ii=1,npel(it1,it2)
625 theta1=zz(1,ipoint)*conv
626 theta2=zz(2,ipoint)*conv
627 ddist=zz(3,ipoint)**3
628 thetap1=zz(4,ipoint)*conv
629 thetap2=zz(5,ipoint)*conv
630 phi=zz(6,ipoint)*conv
631 gam1=zz(7,ipoint)*conv
632 gam2=zz(8,ipoint)*conv
634 print *,"ipoint",ipoint," theta1",theta1," theta2",theta2,
635 & " dd",dd," thetap1",thetap1," thetap2", thetap2," phi",phi,
636 & " gam1",gam1," gam2",gam2
639 facpp(k)=facp1(k)/ddist
642 print *,"facp1",facp1," facpp",facpp
644 call eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
645 & phi,b2m1,b2m,b1m,b1m2,facpp,eneelloc3,eneelloc3b,lder)
647 print *,"eneelloc3",eneelloc3
650 yc(k,ipoint)=eneelloc3(k)
651 diff(k,ipoint)=y(k,ipoint)-yc(k,ipoint)
652 sum=sum+diff(k,ipoint)*diff(k,ipoint)*w(k,ipoint)
655 write (2,'(i7,7f7.1,4(2f10.3,5x))') i,(zz(j,ipoint),j=1,7),
656 & (y(k,ipoint),yc(k,ipoint),k=1,4)
660 c print *,"ipoint",ipoint," l",l
664 if (mask_ello_PMF.gt.0)
665 & jac(mask_ello_PMF)=jac(mask_ello_PMF)+eneelloc3(l)/wello_PMF
669 c jac(2)=eneelloc3(l)/wello
670 c Derivatives in b1 of the first residue
673 jj=mask_bnew1(j,1,it1)
675 if (jj.gt.0) jac(jj)=eneelloc3b(j,1,1,1,l)
676 jj=mask_bnew1(j,2,it1)
679 jac(jj)=eneelloc3b(j,2,1,1,l)
680 else if (it1.lt.0) then
681 jac(jj)=-eneelloc3b(j,2,1,1,l)
685 c print *,"bi1: jac",jac(:n)
686 c Derivatives in b2 of the first residue
689 jj=mask_bnew2(j,1,it1)
690 if (jj.gt.0) jac(jj)=eneelloc3b(j,1,2,1,l)
691 jj=mask_bnew2(j,2,it1)
694 jac(jj)=eneelloc3b(j,2,2,1,l)
695 else if (it1.lt.0) then
696 jac(jj)=-eneelloc3b(j,2,2,1,l)
700 c print *,"bi2: jac",jac(:n)
701 c Derivatives in b1 of the second residue
704 jj=mask_bnew1(j,1,iabs(it2))
705 if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,1,2,l)
706 jj=mask_bnew1(j,2,iabs(it2))
709 jac(jj)=jac(jj)+eneelloc3b(j,2,1,2,l)
710 else if (it2.lt.0) then
711 jac(jj)=jac(jj)-eneelloc3b(j,2,1,2,l)
715 c print *,"bj1: jac",jac(:n)
716 c Derivatives in b2 of the second residue
718 jj=mask_bnew2(j,1,iabs(it2))
719 if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,2,2,l)
720 jj=mask_bnew2(j,2,iabs(it2))
723 jac(jj)=jac(jj)+eneelloc3b(j,2,2,2,l)
724 else if (it2.lt.0) then
725 jac(jj)=jac(jj)-eneelloc3b(j,2,2,2,l)
729 c print *,"bj2: jac",jac(:n)
731 g(k)=g(k)+diff(l,ipoint)*jac(k)*w(l,ipoint)
732 c h(k,k)=h(k,k)+jac(k)**2*w(l,ipoint)
734 c h(k,ll)=h(k,ll)+jac(k)*jac(ll)*w(l,ipoint)
738 print *,"ipoint",ipoint," jac",jac(:n)," w, diff",w(l,ipoint),
739 & diff(1,ipoint)," g",g(:n)
742 c Contributions to the gradient and the hessian
753 c if (nfun.ne.-1) then
755 c write (2,*) "sum",sum
757 call MPI_Reduce(sum_,sum,1,MPI_DOUBLE_PRECISION,
758 & MPI_SUM,Master,Comm1,ierror)
760 c write (iout,*) "after reducing function lder",lder
764 call MPI_Reduce(g_(1),g(1),n,MPI_DOUBLE_PRECISION,
765 & MPI_SUM,Master,Comm1,ierror)
766 c write (iout,*) "after reducing gradient"
771 c if (nfun.ne.-1) nfun=nfun+1
774 c write (iout,*) "Exit rdif"
779 c----------------------------------------------------------------------
780 subroutine eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
781 & phi,bi1,bi2,bj1,bj2,facpp,eneelloc3,eneelloc3b,lder)
784 include "DIMENSIONS.ZSCOPT"
785 include "COMMON.IOUNITS"
786 double precision theta1,theta2,gam1,gam2,dd,thetap1,thetap2,phi
787 double precision bi1(3,2),bi2(3,2),bj1(3,2),bj2(3,2),elpp(2,2),
789 double precision eneelloc3(4),eneelloc3b(3,2,2,2,4),emat1(3,3),
790 & emat2(3,3),emati1(3,3),emati2(3,3),ematj1(3,3),ematj2(3,3)
791 double precision mui1(3),mui2(3),muj1(3),muj2(3)
792 double precision cost1,cost2,sint1,sint2,costp1,sintp1
793 double precision facpp(4)
795 common /calcdip/ cost1,cost2,sint1,sint2
796 double precision pi /3.14159265358979323844d0/
808 print '(4(a,2f10.5))',"bi1",(bi1(j,k),k=1,2)," bi2",
810 & " bj1",(bj1(j,k),k=1,2)," bj2",(bj2(j,k),k=1,2)
812 print *,"cost1",cost1," sint1",sint1," cost2",cost2," sint2",sint2
815 mui1(j+1)=(bi1(1,j)+bi1(2,j)*cost1+bi1(3,j)*cost1*cost1)*sint1
816 muj1(j+1)=(bj1(1,j)+bj1(2,j)*cost2+bj1(3,j)*cost2*cost2)*sint2
817 mui2(j+1)=(bi2(1,j)+bi2(2,j)*cost1+bi2(3,j)*cost1*cost1)*sint1
818 muj2(j+1)=(bj2(1,j)+bj2(2,j)*cost2+bj2(3,j)*cost2*cost2)*sint2
825 print *,"Before rotation"
826 print '(a,3f10.5)',"mui1",mui1
827 print '(a,3f10.5)',"muj1",muj1
828 print '(a,3f10.5)',"mui2",mui2
829 print '(a,3f10.5)',"muj2",muj2
833 c print *,"costp1",costp1," sintp1",sintp1
843 call coorsys(thetap2,phi,emat2)
846 print '(3f10.5)',((emat1(j,k),k=1,3),j=1,3)
848 print '(3f10.5)',((emat2(j,k),k=1,3),j=1,3)
850 call rotmatsub(emat1,gam1,emati1)
851 call rotmatsub(emat1,gam1,emati2)
854 print '(3f10.5)',((emati1(k,j),k=1,3),j=1,3)
856 print '(3f10.5)',((emati2(k,j),k=1,3),j=1,3)
858 call matvecsub(emati1,mui1)
859 call matvecsub(emati2,mui2)
860 call rotmatsub(emat2,gam2,ematj1)
861 call rotmatsub(emat2,gam2,ematj2)
864 print '(3f10.5)',((ematj1(k,j),k=1,3),j=1,3)
866 print '(3f10.5)',((ematj2(k,j),k=1,3),j=1,3)
868 call matvecsub(ematj1,muj1)
869 call matvecsub(ematj2,muj2)
871 print *,"After rotation"
872 print '(a,3f10.5)',"mui1",mui1
873 print '(a,3f10.5)',"muj1",muj1
874 print '(a,3f10.5)',"mui2",mui2
875 print '(a,3f10.5)',"muj2",muj2
877 call edipole(mui1,muj1,emati1,ematj1,facpp(1),eneelloc3(1),
878 & eneelloc3b(1,1,1,1,1),eneelloc3b(1,1,1,2,1),-1,-1,lder)
879 call edipole(mui2,muj1,emati2,ematj1,facpp(2),eneelloc3(2),
880 & eneelloc3b(1,1,2,1,2),eneelloc3b(1,1,1,2,2),1,-1,lder)
881 call edipole(mui1,muj2,emati1,ematj2,facpp(3),eneelloc3(3),
882 & eneelloc3b(1,1,1,1,3),eneelloc3b(1,1,2,2,3),-1,1,lder)
883 call edipole(mui2,muj2,emati2,ematj2,facpp(4),eneelloc3(4),
884 & eneelloc3b(1,1,2,1,4),eneelloc3b(1,1,2,2,4),1,1,lder)
887 c-----------------------------------------------------------------------------
888 subroutine eturn3PMF(theta1,theta2,gam,b1m,b2m,em1,em2,e0m1,e0m2,
889 & facturn,d1,d2,d3,eello_turn3,eello_turn3_b1,eello_turn3_b2,
890 & eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
892 double precision theta1,theta2,gam,cost1,cost2,sint1,sint2,cosg,
893 & sing,sint1sq,sint2sq,aux
894 double precision b1m(3,2),b2m(3,2),b1(2),b2(2),em1(2,2,2),
895 & em2(2,2,2),e0m1(3),e0m2(3),e1(2,2),e2(2,2),facturn,emat1(3,3),
896 & emat2(3,3),emat1t(3,3),mu1(2),mu2(2),e2gam(2,2),eello_turn3,
897 & eello_turn3_b2(3,2),eello_turn3_b1(3,2),eello_turn3_e1(2,2,2),
898 & eello_turn3_e2(2,2,2),eello_turn3_e0m1(3),eello_turn3_e0m2(3),
899 & r1(3),r2(3),er(3),a(3,3),atemp(3,3),etemp1(2,2),etemp2(2,2),
901 double precision pi /3.14159265358979323844d0/
916 mu1(1)=(b1m(1,1)+b1m(2,1)*cost1+b1m(3,1)*cost1*cost1)*sint1
917 mu1(2)=(b1m(1,2)+b1m(2,2)*cost1+b1m(3,2)*cost1*cost1)*sint1
918 mu2(1)=(b2m(1,1)+b2m(2,1)*cost2+b2m(3,1)*cost2*cost2)*sint2
919 mu2(2)=(b2m(1,2)+b2m(2,2)*cost2+b2m(3,2)*cost2*cost2)*sint2
924 e1(k,l)=(em1(1,k,l)+em1(2,k,l)*cost1)*sint1sq
925 e2(k,l)=(em2(1,k,l)+em2(2,k,l)*cost2)*sint2sq
928 e1(1,1)=e1(1,1)+e0m1(1)*cost1
929 e1(1,2)=e1(1,2)+e0m1(2)+e0m1(3)*cost1
930 e1(2,1)=e1(2,1)+e0m1(2)*cost1+e0m1(3)
931 e1(2,2)=e1(2,2)-e0m1(1)
932 e2(1,1)=e2(1,1)+e0m2(1)*cost2
933 e2(1,2)=e2(1,2)+e0m2(2)+e0m2(3)*cost2
934 e2(2,1)=e2(2,1)+e0m2(2)*cost2+e0m2(3)
935 e2(2,2)=e2(2,2)-e0m2(1)
937 print *,"mu1",mu1," mu2",mu2," e1",e1," e2",e2
939 c Coordinate system of the first dipole
951 emat1t(i,j)=emat1(j,i)
955 print *,"theta1",theta1*180/pi," theta2",theta2*180/pi,
958 c Calculate the reference system of the second unit
959 call coorsys1(theta2,gam,emat2)
960 c Dipole-dipole interaction matrix
962 print *,"d1",d1," d2",d2," d3",d3
963 print *,"emat1",emat1," emat2",emat2
968 call matvecsub(emat1,r1)
975 call matvecsub(emat2,r2)
980 ddist=dsqrt((r2(1)-r1(1))**2+(r2(2)-r1(2))**2+(r2(3)-r1(3))**2)
982 er(1)=(r2(1)-r1(1))/ddist
983 er(2)=(r2(2)-r1(2))/ddist
984 er(3)=(r2(3)-r1(3))/ddist
986 print *,"norm",er(1)**2+er(2)**2+er(3)**2
987 print *,"r1",r1," r2",r2
989 print *," dd",dd," dd3",dd3
993 a(i,j)=-3*er(i)*er(j)*dd3
994 c a(i,j)=-3*er(i)*er(j)
999 c a(i,i)=1.0d0+a(i,i)
1002 print *,"a before transformation",a
1003 print *,"emat1t",emat1t
1004 print *,"emat2",emat2
1006 call matmat(emat1t,a,atemp)
1008 print *,"a after transformation",atemp
1010 call matmat(atemp,emat2,a)
1012 print *,"a after transformation",a
1014 print "(1h[,2f10.5,1h],5x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
1015 & mu1(1),mu1(2),a(2,2),a(2,3),mu2(1)
1016 print "(27x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
1017 & a(3,2),a(3,3),mu2(2)
1019 eello_turn3=(a(2,2)*mu1(1)*mu2(1)+a(2,3)*mu1(1)*mu2(2)
1020 & +a(3,2)*mu1(2)*mu2(1)+a(3,3)*mu1(2)*mu2(2))
1022 print *,"eello_turn3_1",eello_turn3
1024 eello_turn3_b1(1,1)=-(a(2,2)*mu2(1)+a(2,3)*mu2(2))*sint1
1025 eello_turn3_b1(2,1)=eello_turn3_b1(1,1)*cost1
1026 eello_turn3_b1(3,1)=eello_turn3_b1(2,1)*cost1
1027 eello_turn3_b1(1,2)=(a(3,2)*mu2(1)+a(3,3)*mu2(2))*sint1
1028 eello_turn3_b1(2,2)=eello_turn3_b1(1,2)*cost1
1029 eello_turn3_b1(3,2)=eello_turn3_b1(2,2)*cost1
1030 eello_turn3_b2(1,1)=-(a(2,2)*mu1(1)+a(3,2)*mu1(2))*sint2
1031 eello_turn3_b2(2,1)=eello_turn3_b2(1,1)*cost2
1032 eello_turn3_b2(3,1)=eello_turn3_b2(2,1)*cost2
1033 eello_turn3_b2(1,2)=(a(2,3)*mu1(1)+a(3,3)*mu1(2))*sint2
1034 eello_turn3_b2(2,2)=eello_turn3_b2(1,2)*cost2
1035 eello_turn3_b2(3,2)=eello_turn3_b2(2,2)*cost2
1036 e2gam(1,1)=e2(1,1)*cosg+e2(2,1)*sing
1037 e2gam(1,2)=e2(1,2)*cosg+e2(2,2)*sing
1038 e2gam(2,1)=-e2(1,1)*sing+e2(2,1)*cosg
1039 e2gam(2,2)=-e2(1,2)*sing+e2(2,2)*cosg
1054 aux=aux+a(i+1,j+1)*(e1(i,1)*e2(1,j)+e1(i,2)*e2(2,j))
1058 print *,"eello_turn3_2",aux/2
1060 eello_turn3=eello_turn3+aux/2
1062 c Derivatives in the elements of complete e1 and e2 matrices
1065 etemp1(i,j)=(a(i+1,2)*e2(j,1)+a(i+1,3)*e2(j,2))/2
1066 etemp2(i,j)=(a(2,j+1)*e1(1,i)+a(3,j+1)*e1(2,i))/2
1069 c Derivatives in the components of e1
1071 eello_turn3_e1(1,1,i)=etemp1(1,i)*sint1sq
1072 eello_turn3_e1(2,1,i)=eello_turn3_e1(1,1,i)*cost1
1073 eello_turn3_e1(1,2,i)= etemp1(2,i)*sint1sq
1074 eello_turn3_e1(2,2,i)=eello_turn3_e1(1,2,i)*cost1
1076 c Derivatives in the components of e01m
1077 eello_turn3_e0m1(1)=etemp1(1,1)*cost1-etemp1(2,2)
1078 eello_turn3_e0m1(2)=etemp1(1,2)+cost1*etemp1(2,1)
1079 eello_turn3_e0m1(3)=etemp1(1,2)*cost1+etemp1(2,1)
1080 c Derivatives in the elements of e2
1082 eello_turn3_e2(1,1,i)=-(etemp2(1,i)*cosg+etemp2(2,i)*sing)
1084 eello_turn3_e2(2,1,i)=eello_turn3_e2(1,1,i)*cost2
1085 eello_turn3_e2(1,2,i)=(-etemp2(1,i)*sing+etemp2(2,i)*cosg)
1087 eello_turn3_e2(2,2,i)=eello_turn3_e2(1,2,i)*cost2
1089 c Derivatives in the elements of e02m
1090 eello_turn3_e0m2(1)=-(etemp2(1,1)*cosg+etemp2(2,1)*sing)*cost2
1091 & +etemp2(1,2)*sing-etemp2(2,2)*cosg
1092 eello_turn3_e0m2(2)=-etemp2(1,2)*cosg-etemp2(2,2)*sing+
1093 & (-etemp2(1,1)*sing+etemp2(2,1)*cosg)*cost2
1094 eello_turn3_e0m2(3)=(-etemp2(1,2)*cosg-etemp2(2,2)*sing)*cost2
1095 & -etemp2(1,1)*sing+etemp2(2,1)*cosg
1098 c-----------------------------------------------------------------------------
1099 subroutine edipole(mui,muj,emati,ematj,facpp,ene,enebi,enebj,
1101 include "DIMENSIONS"
1102 include "DIMENSIONS.ZSCOPT"
1103 include "COMMON.IOUNITS"
1104 double precision mui(3),muj(3),emati(3,3),ematj(3,3),
1105 & facpp,ene,enebi(3,2),enebj(3,2),muimuj,cost1,sint1,cost2,sint2
1108 common /calcdip/ cost1,cost2,sint1,sint2
1109 muimuj=mui(1)*muj(1)+mui(2)*muj(2)+mui(3)*muj(3)
1110 ene=facpp*(muimuj-3*mui(3)*muj(3))
1111 if (idiri.gt.0) then
1112 enebi(1,1)=(emati(1,2)*muj(1)+emati(2,2)*muj(2)
1113 & +emati(3,2)*muj(3)-3*muj(3)*emati(3,2))*sint1*facpp
1114 enebi(2,1)=enebi(1,1)*cost1
1115 enebi(3,1)=enebi(2,1)*cost1
1116 enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
1117 & +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
1118 enebi(2,2)=enebi(1,2)*cost1
1119 enebi(3,2)=enebi(2,2)*cost1
1121 enebi(1,1)=(-emati(1,2)*muj(1)-emati(2,2)*muj(2)
1122 & -emati(3,2)*muj(3)+3*muj(3)*emati(3,2))*sint1*facpp
1123 enebi(2,1)=enebi(1,1)*cost1
1124 enebi(3,1)=enebi(2,1)*cost1
1125 enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
1126 & +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
1127 enebi(2,2)=enebi(1,2)*cost1
1128 enebi(3,2)=enebi(2,2)*cost1
1130 if (idirj.gt.0) then
1131 enebj(1,1)=(ematj(1,2)*mui(1)+ematj(2,2)*mui(2)
1132 & +ematj(3,2)*mui(3)-3*mui(3)*ematj(3,2))*sint2*facpp
1133 enebj(2,1)=enebj(1,1)*cost2
1134 enebj(3,1)=enebj(2,1)*cost2
1135 enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
1136 & +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
1137 enebj(2,2)=enebj(1,2)*cost2
1138 enebj(3,2)=enebj(2,2)*cost2
1140 enebj(1,1)=(-ematj(1,2)*mui(1)-ematj(2,2)*mui(2)
1141 & -ematj(3,2)*mui(3)+3*mui(3)*ematj(3,2))*sint2*facpp
1142 enebj(2,1)=enebj(1,1)*cost2
1143 enebj(3,1)=enebj(2,1)*cost2
1144 enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
1145 & +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
1146 enebj(2,2)=enebj(1,2)*cost2
1147 enebj(3,2)=enebj(2,2)*cost2
1151 c-----------------------------------------------------------------------------
1152 subroutine matvecsub(mat,vec)
1154 double precision mat(3,3),vec(3),auxvec(3)
1159 auxvec(i)=auxvec(i)+mat(i,j)*vec(j)
1165 c-------------------------------------------------------------------------------
1166 subroutine rotmatsub(mat,gam,matgam)
1168 double precision mat(3,3),gam,matgam(3,3),cosg,sing
1173 matgam(i,1)=mat(i,1)
1174 matgam(i,2)= mat(i,2)*cosg+mat(i,3)*sing
1175 matgam(i,3)=-mat(i,2)*sing+mat(i,3)*cosg
1179 c-----------------------------------------------------------------------------
1180 subroutine coorsys(theta,phi,emat)
1182 include "DIMENSIONS"
1183 include "DIMENSIONS.ZSCOPT"
1184 include "COMMON.IOUNITS"
1185 double precision theta,phi,emat(3,3)
1186 double precision cost,sint,cosphi,sinphi
1187 c print *,"theta",theta," phi",phi
1193 print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1196 emat(1,1)=sint*cosphi
1197 emat(2,1)=sint*sinphi
1202 emat(1,3)=-cost*cosphi
1203 emat(2,3)=-cost*sinphi
1207 c-----------------------------------------------------------------------------
1208 subroutine coorsys1(theta,phi,emat)
1210 include "DIMENSIONS"
1211 include "DIMENSIONS.ZSCOPT"
1212 include "COMMON.IOUNITS"
1213 double precision theta,phi,emat(3,3)
1214 double precision cost,sint,cosphi,sinphi
1215 c print *,"theta",theta," phi",phi
1221 print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1224 emat(1,1)=sint*cosphi
1226 emat(3,1)=-sint*sinphi
1227 emat(1,2)=-cost*cosphi
1229 emat(3,2)=cost*sinphi
1235 c-------------------------------------------------------------------------------
1236 subroutine matmat(a,b,c)
1238 double precision a(3,3),b(3,3),c(3,3)
1244 c(i,j)=c(i,j)+a(i,k)*b(k,j)