source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/funder_eello.F

   1       double precision function rdif(n,x,g,ider)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       integer n,nmax
   6       parameter (nmax=max_paropt)
   7 #ifdef MPI
   8       include 'mpif.h'
   9       include 'COMMON.MPI'
  10       integer ierror
  11       double precision g_(max_paropt),sum_
  12 #endif
  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),
  19      &  eello_turn3,
  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
  35       logical lder
  36       double precision x(n),g(n)
  37       integer ipoint,iipoint,ip,it1,it2,itp1,itp2,i,j,k,l,ii,jj!,irt1,irt2
  38       integer ider
  39 #ifdef MPI
  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
  43 c of computations
  44 c      if (nfun.ne.-1) then
  45 c      call MPI_Bcast(n,1,MPI_INTEGER,Master,MPI_COMM_WORLD,ierror)
  46 c      write (2,*) "n",n
  47 c      call flush(iout)
  48       call MPI_Bcast(ider,1,MPI_LOGICAL,Master,MPI_COMM_WORLD,ierror)
  49       if (ider.le.0) return
  50 c      call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION,Master,
  51 c     &    MPI_COMM_WORLD,ierror)
  52 c      endif
  53 c      write (2,*) "x",(x(i),i=1,n)
  54 #endif
  55       lder=ider.gt.1
  56 c      write (iout,*) "ider",ider," lder",lder
  57 c      call flush(iout)
  58 c      print *,"rdif npoint",npoint
  59       IF (lder) THEN
  60 c Zero out gradient and hessian
  61       do i=1,n
  62         g(i)=0.0d0
  63       enddo
  64       ENDIF
  65
  66       ipoint=0
  67       ip=0
  68       b=0.0d0
  69       sum=0.0d0
  70
  71 c      x(16)=dabs(x(16))
  72 c      wello=x(16)
  73 c-- Following for ala-ala only
  74 c      x(2)=dabs(x(2))
  75 c      wello=x(2)
  76 c      wturn=x(3)
  77 c      x(16)=dabs(x(16))
  78 c      x(17)=dabs(x(17))
  79 c      wello=x(16)
  80 c#ifndef TORUN
  81 c      wello=1.0d0
  82 c#endif
  83 c      wturn=x(17)
  84
  85       do it1=0,nloctyp-1
  86       do it2=-nloctyp+1,nloctyp-1
  87 #ifdef DEBUG
  88       write (iout,*) "it1",it1," it2",it2," e0",e0(it1,it2)
  89 #endif
  90 c following to test with selected types
  91 c      do it1=1,1
  92 c      do it2=-1,-1
  93
  94       if (it1.lt.nloctyp-1) then
  95         itp1=1
  96       else
  97         itp1=2
  98       endif
  99       if (iabs(it2).lt.nloctyp-1) then
 100         itp2=1
 101       else
 102         itp2=2
 103       endif
 104       d1=dCACA(itp1)
 105       d2=dCACA(itp2)
 106       d3=dCACA(1)
 107       if (itp1.eq.1) then
 108         facturn3=wturn_PMF*dsqrt(epspp(1))*4.2d0**3
 109       else
 110         facturn3=wturn_PMF*dsqrt(epspp(2))*4.2d0**3
 111       endif
 112
 113       ip=ip+1
 114
 115       delta=180.0d0
 116 c      e0=x(ip)
 117
 118       facp1(4)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
 119       if (itp2.eq.1) then
 120         facp1(2)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
 121       else
 122         facp1(2)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
 123       endif
 124       if (itp1.eq.1) then
 125         facp1(3)=wello_PMF*dsqrt(epspp(1))*4.2d0**3
 126       else
 127         facp1(3)=wello_PMF*dsqrt(epspp(2))*4.2d0**3
 128       endif
 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
 133       else
 134         facp1(1)=wello_PMF*dsqrt(epspp(3))*4.2d0**3
 135       endif
 136 c      print *,"itp1",itp1," itp2",itp2," facp1",facp1
 137
 138 c      irt1=15+20*it1
 139 c      irt2=15+20*iabs(it2)
 140 c-- The following for ALL residues
 141 c      irt1=17+31*it1
 142 c      irt2=17+31*iabs(it2)
 143 c-- The followinig for ala-ala only
 144 c      irt1=3
 145 c      irt2=3
 146 c       print *,"it1",it1," it2",it2
 147 c       print *,irt1,irt2,irt1+31,irt2+31
 148 c----------------------
 149 c      if (FOURIERSPLIT) then
 150         do i=1,2
 151           do k=1,3
 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)
 156           enddo
 157         enddo
 158         do i=1,2
 159           do k=1,2
 160             cm(k,i)=ccnewtor(k,i,it2)
 161             dm(k,i)=ddnewtor(k,i,it1)
 162           enddo
 163         enddo
 164         do i=1,2
 165           do j=1,2
 166             do k=1,2
 167               em(k,j,i)=eenewtor(k,j,i,it1)
 168               em2(k,j,i)=eenewtor(k,j,i,it2)
 169             enddo
 170           enddo
 171         enddo
 172         do k=1,3
 173           e0m(k)=e0newtor(k,it1)
 174           e0m2(k)=e0newtor(k,it2)
 175         enddo
 176 c Invert coefficients for L-chirality
 177 #ifdef DEBUG
 178       write(iout,*) "it1",it1," it2",it2
 179       write(iout,*) "b1i (b1m)"
 180       do i=1,3
 181         write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
 182       enddo
 183       write(iout,*) "b2j (b1m2)"
 184       do i=1,3
 185         write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
 186       enddo
 187       write(iout,*) "b1j (b2m1)"
 188       do i=1,3
 189         write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
 190       enddo
 191       write(iout,*) "b2j (b2m)"
 192       do i=1,3
 193         write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
 194       enddo
 195       write (iout,*) "c"
 196       do i=1,2
 197         write (iout,'(2f10.5)') (cm(i,j),j=1,2)
 198       enddo
 199       write (iout,*) "d"
 200       do i=1,2
 201         write (iout,'(2f10.5)') (dm(i,j),j=1,2)
 202       enddo
 203       print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
 204 #endif
 205       v1_kcc_loc=0.0d0
 206       v2_kcc_loc=0.0d0
 207       etoriv1=0.0d0
 208       etoriv2=0.0d0
 209 c      cm=0.0d0
 210 c      dm=0.0d0
 211       do k=1,3
 212         do l=1,3
 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)
 215         enddo
 216       enddo
 217       do k=1,2
 218         do l=1,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))
 221         enddo
 222       enddo
 223 #define CHUJUN
 224 #ifdef CHUJUN
 225       if (torsion_pmf) then
 226
 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)
 235       iipoint=ipoint
 236       do ii=1,np(it1,it2)
 237         ipoint=ipoint+1
 238         i=ipoint
 239         theta1=zz(1,i)*conv
 240         theta2=zz(2,i)*conv
 241         gam=(zz(3,i)+delta-180)*conv
 242         cost1=dcos(theta1)
 243         cost2=dcos(theta2)
 244         sint1=dsin(theta1)
 245         sint2=dsin(theta2)
 246         cosg=dcos(gam)
 247         sing=dsin(gam)
 248         cos2g=dcos(2*gam)
 249         sin2g=dsin(2*gam)
 250 #ifdef OLDCALC
 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
 259 #ifdef DEBUG
 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,
 263      &    " sin2g",sin2g
 264 #endif
 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
 275 #endif
 276 c check with the formula from unres
 277         c1(0)=0.0d0
 278         c2(0)=0.0d0
 279         c1(1)=1.0d0
 280         c2(1)=1.0d0
 281         do j=2,3
 282           c1(j)=c1(j-1)*cost1
 283           c2(j)=c2(j-1)*cost2
 284         enddo
 285         etori=e0(it1,it2)
 286         sint1sint2=1.0d0
 287         do j=1,2
 288           cosphi=dcos(j*gam)
 289           sinphi=dsin(j*gam)
 290           sint1sint2=sint1sint2*sint1*sint2
 291           sumvalc=0.0d0
 292           sumvals=0.0d0
 293           do k=1,3
 294             do l=1,3
 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
 299             enddo
 300           enddo
 301 c          if (j.eq.1) print *,"j",j," sumvlac",sumvalc*sint1sint2,
 302 c     &      b21*b11+b22*b12,
 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
 308         enddo
 309 c      do k=1,3
 310 c        do l=1,3
 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)
 313 c        enddo
 314 c      enddo
 315 c      do k=1,2
 316 c        do l=1,2
 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))
 319 c        enddo
 320 c      enddo
 321 c        print *,"i",i," theta",zz(i,1),zz(i,2)," gamma",zz(i,3),
 322 c     &    " f",f," etori",etori
 323 #define CIPUN
 324 #ifdef CIPUN
 325         yc(1,i)=etori
 326 #ifdef DEBUG
 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)
 329 #endif
 330         diff(1,i)=y(1,i)-yc(1,i)
 331         sum=sum+diff(1,i)*diff(1,i)*w(1,i)
 332 c
 333 c Compute derivatives of the energy in torsional parameters
 334 c
 335         IF (lder) THEN
 336 c
 337         do k=1,n
 338           jac(k)=0.0d0
 339         enddo
 340 c Free term
 341         if (mask_e0(it1,it2).gt.0) jac(mask_e0(it1,it2))=1.0d0
 342 c Derivatives in b1 of the second residue
 343 c        jj=-1
 344         do j=1,3
 345 c          jj=jj+2
 346           jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
 347           if (jj.gt.0) then
 348             do k=1,3
 349               jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
 350      &        +etoriv2(j,k,1)*b2m(k,2)
 351             enddo
 352           endif
 353           jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
 354           if (jj.gt.0) then
 355             do k=1,3
 356               if (it2.gt.0) then
 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)
 362               endif
 363             enddo
 364           endif
 365         enddo
 366 c Derivatives in b2 of the first residue
 367         do j=1,3
 368 c          jj=jj+2
 369           jj=iabs(mask_bnew2tor(j,1,it1))
 370           if (jj.gt.0) then
 371             do k=1,3
 372               jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
 373      &        +etoriv2(k,j,1)*b1m(k,2)
 374             enddo
 375           endif
 376           jj=iabs(mask_bnew2tor(j,2,it1))
 377           if (jj.gt.0) then
 378             do k=1,3
 379               if (it1.gt.0) then
 380                 jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
 381      &          +etoriv2(k,j,1)*b1m(k,1)
 382               endif
 383             enddo
 384           endif
 385         enddo
 386 c Derivatives in c
 387         do j=1,2
 388 c          jj=jj+2
 389           jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
 390           if (jj.gt.0) then
 391             do k=1,2
 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)
 396             enddo
 397           endif
 398           jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
 399           if (jj.gt.0) then
 400             do k=1,2
 401               if (it2.gt.0) then
 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)
 411               endif
 412             enddo
 413           endif
 414         enddo
 415 c Derivatives in d
 416         do j=1,2
 417 c          jj=jj+2
 418           jj=iabs(mask_ddnewtor(j,1,it1))
 419           if (jj.gt.0) then
 420             do k=1,2
 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)
 425             enddo
 426           endif
 427           jj=iabs(mask_ddnewtor(j,2,it1))
 428           if (jj.gt.0) then
 429             do k=1,2
 430               if (it1.gt.0) then
 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)
 435               endif
 436             enddo
 437           endif
 438         enddo
 439 c Contributions to the gradient and the hessian
 440         do k=1,n
 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)
 443 c          do l=1,k-1
 444 c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(1,i)
 445 c          enddo
 446         enddo
 447 #ifdef DEBUG
 448         print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
 449      &   " g",g(:n)
 450 #endif
 451         ENDIF ! lder
 452 #endif
 453       enddo  ! ii
 454
 455       endif ! torsion_pmf
 456
 457       do i=1,2
 458         do k=1,3
 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)
 463         enddo
 464       enddo
 465       do i=1,2
 466         do j=1,2
 467           do k=1,2
 468             em(k,j,i)=eenew(k,j,i,it1)
 469             em2(k,j,i)=eenew(k,j,i,it2)
 470           enddo
 471         enddo
 472       enddo
 473       do k=1,3
 474         e0m(k)=e0new(k,it1)
 475         e0m2(k)=e0new(k,it2)
 476       enddo
 477
 478       if (turn_pmf) then
 479
 480       ipoint=iipoint
 481       do ii=1,np(it1,it2)
 482         ipoint=ipoint+1
 483         i=ipoint
 484         theta1=zz(1,i)*conv
 485         theta2=zz(2,i)*conv
 486         gam=(zz(3,i)+delta-180)*conv
 487 #define CYCUN
 488 #ifdef CYCUN
 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)
 493         yc(2,i)=eello_turn3
 494 #ifdef DEBUG
 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)
 497 #endif
 498         diff(2,i)=y(2,i)-yc(2,i)
 499         sum=sum+diff(2,i)*diff(2,i)*w(2,i)
 500
 501         IF (lder) THEN
 502
 503         jac=0.0d0
 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
 506 c        jj=-1
 507         do j=1,3
 508           jj=mask_bnew1(j,1,it1)
 509           if (jj.gt.0) then
 510 c            jj=jj+2
 511             jac(jj)=jac(jj)+eello_turn3_b2(j,1)
 512           endif
 513         enddo
 514         do j=1,3
 515           jj=mask_bnew1(j,2,it1)
 516           if (jj.gt.0) then
 517             if (it1.ne.0) then
 518               jac(jj)=jac(jj)+eello_turn3_b2(j,2)
 519             endif
 520           endif
 521         enddo
 522 c Derivatives of in b2 of the second residue
 523         do j=1,3
 524           jj=mask_bnew2(j,1,iabs(it2))
 525           if (jj.gt.0) then
 526             jac(jj)=jac(jj)+eello_turn3_b1(j,1)
 527           endif
 528         enddo
 529         do j=1,3
 530           jj=mask_bnew2(j,2,iabs(it2))
 531           if (jj.gt.0) then
 532             if (it2.gt.0) then
 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)
 536             endif
 537           endif
 538         enddo
 539 #define EE
 540 #ifdef EE
 541 c Derivatives in e of the first residue
 542 c        jj=21
 543         do j=1,2
 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)
 548           if (it1.ne.0) then
 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)
 553           endif
 554 c          jj=jj+4
 555         enddo
 556         jj=mask_e0new(1,it1)
 557         if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(1)
 558         if (it1.ne.0) then
 559           jj=mask_e0new(2,it1)
 560           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(2)
 561           jj=mask_e0new(3,it1)
 562           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(3)
 563         endif
 564 c Derivatives in e of the second residue
 565 c        jj=21
 566         do j=1,2
 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)
 571           if (it2.gt.0) then
 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)
 581           endif
 582 c          jj=jj+4
 583         enddo
 584         jj=mask_e0new(1,iabs(it2))
 585         if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(1)
 586         if (it2.gt.0) then
 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)
 596         endif
 597 #endif
 598 c Contributions to the gradient and the hessian
 599         do k=1,n
 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)
 602 c          do l=1,k-1
 603 c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(2,i)
 604 c          enddo
 605         enddo
 606 #ifdef DEBUG
 607         print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
 608      &   " g",g(:n)
 609 #endif
 610 c
 611         ENDIF ! lder
 612 #endif
 613       enddo
 614
 615       endif ! turn_pmf
 616 c      nfun=nfun+1
 617 c      rdif=sum
 618 c      return
 619 #endif
 620       if (eello_pmf) then
 621
 622       do ii=1,npel(it1,it2)
 623         eneelloc3b=0.0d0
 624         ipoint=ipoint+1
 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
 633 #ifdef DEBUG
 634         print *,"ipoint",ipoint," theta1",theta1," theta2",theta2,
 635      &   " dd",dd," thetap1",thetap1," thetap2", thetap2," phi",phi,
 636      &   " gam1",gam1," gam2",gam2
 637 #endif
 638         do k=1,4
 639           facpp(k)=facp1(k)/ddist
 640         enddo
 641 #ifdef DEBUG
 642         print *,"facp1",facp1," facpp",facpp
 643 #endif
 644         call eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
 645      &    phi,b2m1,b2m,b1m,b1m2,facpp,eneelloc3,eneelloc3b,lder)
 646 #ifdef DEBUG
 647         print *,"eneelloc3",eneelloc3
 648 #endif
 649         do k=1,4
 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)
 653         enddo
 654 #ifdef DEBUG
 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)
 657 #endif
 658         IF (lder) THEN
 659         do l=1,4
 660 c          print *,"ipoint",ipoint," l",l
 661           do k=1,n
 662             jac(k)=0.0d0
 663           enddo
 664           if (mask_ello_PMF.gt.0)
 665      &    jac(mask_ello_PMF)=jac(mask_ello_PMF)+eneelloc3(l)/wello_PMF
 666 c#ifndef TORUN
 667 c          jac(16)=0.0d0
 668 c#endif
 669 c          jac(2)=eneelloc3(l)/wello
 670 c Derivatives in b1 of the first residue
 671 c          jj=-1
 672           do j=1,3
 673             jj=mask_bnew1(j,1,it1)
 674 c            jj=jj+2
 675             if (jj.gt.0) jac(jj)=eneelloc3b(j,1,1,1,l)
 676             jj=mask_bnew1(j,2,it1)
 677             if (jj.gt.0) then
 678               if (it1.gt.0) then
 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)
 682               endif
 683             endif
 684           enddo
 685 c          print *,"bi1: jac",jac(:n)
 686 c Derivatives in b2 of the first residue
 687           do j=1,3
 688 c            jj=jj+2
 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)
 692             if (jj.gt.0) then
 693               if (it1.gt.0) then
 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)
 697               endif
 698             endif
 699           enddo
 700 c          print *,"bi2: jac",jac(:n)
 701 c Derivatives in b1 of the second residue
 702 c          jj=-1
 703           do j=1,3
 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))
 707             if (jj.gt.0) then
 708               if (it2.gt.0) then
 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)
 712               endif
 713             endif
 714           enddo
 715 c          print *,"bj1: jac",jac(:n)
 716 c Derivatives in b2 of the second residue
 717           do j=1,3
 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))
 721             if (jj.gt.0) then
 722               if (it2.gt.0) then
 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)
 726               endif
 727             endif
 728           enddo
 729 c          print *,"bj2: jac",jac(:n)
 730           do k=1,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)
 733 c            do ll=1,k-1
 734 c              h(k,ll)=h(k,ll)+jac(k)*jac(ll)*w(l,ipoint)
 735 c            enddo
 736           enddo
 737 #ifdef DEBUG
 738           print *,"ipoint",ipoint," jac",jac(:n)," w, diff",w(l,ipoint),
 739      &     diff(1,ipoint)," g",g(:n)
 740 #endif
 741         enddo ! l
 742 c Contributions to the gradient and the hessian
 743
 744         ENDIF
 745       enddo ! ii
 746
 747       endif ! eello_PMF
 748
 749       enddo ! it2
 750       enddo ! it1
 751
 752 #ifdef MPI
 753 c      if (nfun.ne.-1) then
 754       sum_=sum
 755 c      write (2,*) "sum",sum
 756       call flush(2)
 757       call MPI_Reduce(sum_,sum,1,MPI_DOUBLE_PRECISION,
 758      &   MPI_SUM,Master,Comm1,ierror)
 759 c      endif
 760 c      write (iout,*) "after reducing function lder",lder
 761 c      call flush(iout)
 762       if (lder) then
 763         g_=g
 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"
 767 c      call flush(iout)
 768       endif
 769 #endif
 770
 771 c      if (nfun.ne.-1) nfun=nfun+1
 772
 773       rdif=0.5d0*sum
 774 c      write (iout,*) "Exit rdif"
 775 c      call flush(iout)
 776 c      stop
 777       return
 778       end
 779 c----------------------------------------------------------------------
 780       subroutine eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
 781      &    phi,bi1,bi2,bj1,bj2,facpp,eneelloc3,eneelloc3b,lder)
 782       implicit none
 783       include "DIMENSIONS"
 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),
 788      &    r0pp(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)
 794       logical lder
 795       common /calcdip/ cost1,cost2,sint1,sint2
 796       double precision pi /3.14159265358979323844d0/
 797       integer j,k
 798       mui1=0.0d0
 799       mui2=0.0d0
 800       muj1=0.0d0
 801       muj2=0.0d0
 802       cost1=dcos(theta1)
 803       sint1=dsin(theta1)
 804       cost2=dcos(theta2)
 805       sint2=dsin(theta2)
 806 #ifdef DEBUG
 807       do j=1,3
 808       print '(4(a,2f10.5))',"bi1",(bi1(j,k),k=1,2)," bi2",
 809      &  (bi2(j,k),k=1,2),
 810      & " bj1",(bj1(j,k),k=1,2)," bj2",(bj2(j,k),k=1,2)
 811       enddo
 812       print *,"cost1",cost1," sint1",sint1," cost2",cost2," sint2",sint2
 813 #endif
 814       do j=1,2
 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
 819       enddo
 820 c      mui2(2)=-mui2(2)
 821 c      muj2(2)=-muj2(2)
 822       mui1(2)=-mui1(2)
 823       muj1(2)=-muj1(2)
 824 #ifdef DEBUG
 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
 830 #endif
 831       sintp1=dsin(thetap1)
 832       costp1=dcos(thetap1)
 833 c      print *,"costp1",costp1," sintp1",sintp1
 834       emat1(1,1)=sintp1
 835       emat1(1,2)=0.0d0
 836       emat1(1,3)=-costp1
 837       emat1(2,1)=0.0d0
 838       emat1(2,2)=1.0d0
 839       emat1(2,3)=0.0d0
 840       emat1(3,1)=costp1
 841       emat1(3,2)=0.0d0
 842       emat1(3,3)=sintp1
 843       call coorsys(thetap2,phi,emat2)
 844 #ifdef DEBUG
 845       print *,"emat1"
 846       print '(3f10.5)',((emat1(j,k),k=1,3),j=1,3)
 847       print *,"emat2"
 848       print '(3f10.5)',((emat2(j,k),k=1,3),j=1,3)
 849 #endif
 850       call rotmatsub(emat1,gam1,emati1)
 851       call rotmatsub(emat1,gam1,emati2)
 852 #ifdef DEBUG
 853       print *,"emati1"
 854       print '(3f10.5)',((emati1(k,j),k=1,3),j=1,3)
 855       print *,"emati2"
 856       print '(3f10.5)',((emati2(k,j),k=1,3),j=1,3)
 857 #endif
 858       call matvecsub(emati1,mui1)
 859       call matvecsub(emati2,mui2)
 860       call rotmatsub(emat2,gam2,ematj1)
 861       call rotmatsub(emat2,gam2,ematj2)
 862 #ifdef DEBUG
 863       print *,"ematj1"
 864       print '(3f10.5)',((ematj1(k,j),k=1,3),j=1,3)
 865       print *,"ematj2"
 866       print '(3f10.5)',((ematj2(k,j),k=1,3),j=1,3)
 867 #endif
 868       call matvecsub(ematj1,muj1)
 869       call matvecsub(ematj2,muj2)
 870 #ifdef DEBUG
 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
 876 #endif
 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)
 885       return
 886       end
 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)
 891       implicit none
 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),
 900      &  d1,d2,d3,ddist,dd3
 901       double precision pi /3.14159265358979323844d0/
 902       integer i,j,k,l
 903 c diagnostics
 904 c      d1=3.8
 905 c      d2=3.8
 906 c      d3=3.8
 907 c end diagnostics
 908       cost1=dcos(theta1)
 909       cost2=dcos(theta2)
 910       sint1=dsin(theta1)
 911       sint1sq=sint1*sint1
 912       sint2=dsin(theta2)
 913       sint2sq=sint2*sint2
 914       cosg=dcos(gam)
 915       sing=dsin(gam)
 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
 920       mu1(1)=-mu1(1)
 921       mu2(1)=-mu2(1)
 922       do k=1,2
 923         do l=1,2
 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
 926         enddo
 927       enddo
 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)
 936 #ifdef DEBUG
 937       print *,"mu1",mu1," mu2",mu2," e1",e1," e2",e2
 938 #endif
 939 c Coordinate system of the first  dipole
 940       emat1(1,1)=sint1
 941       emat1(1,2)=cost1
 942       emat1(1,3)=0.0d0
 943       emat1(2,1)=-cost1
 944       emat1(2,2)=sint1
 945       emat1(2,3)=0.0d0
 946       emat1(3,1)=0.0d0
 947       emat1(3,2)=0.0d0
 948       emat1(3,3)=1.0d0
 949       do i=1,3
 950         do j=1,3
 951           emat1t(i,j)=emat1(j,i)
 952         enddo
 953       enddo
 954 #ifdef DEBUG
 955       print *,"theta1",theta1*180/pi," theta2",theta2*180/pi,
 956      &  " gam",gam*180/pi
 957 #endif
 958 c Calculate the reference system of the second unit
 959       call coorsys1(theta2,gam,emat2)
 960 c Dipole-dipole interaction matrix
 961 #ifdef DEBUG
 962       print *,"d1",d1," d2",d2," d3",d3
 963       print *,"emat1",emat1," emat2",emat2
 964 #endif
 965       r1(1)=-d1/2
 966       r1(2)=0.0d0
 967       r1(3)=0.0d0
 968       call matvecsub(emat1,r1)
 969       r2(1)=-d3/2
 970       r2(2)=0.0d0
 971       r2(3)=0.0d0
 972 #ifdef DEBUG
 973       print *," r2",r2
 974 #endif
 975       call matvecsub(emat2,r2)
 976 #ifdef DEBUG
 977       print *," r2",r2
 978 #endif
 979       r2(2)=r2(2)+d2
 980       ddist=dsqrt((r2(1)-r1(1))**2+(r2(2)-r1(2))**2+(r2(3)-r1(3))**2)
 981       dd3=facturn/ddist**3
 982       er(1)=(r2(1)-r1(1))/ddist
 983       er(2)=(r2(2)-r1(2))/ddist
 984       er(3)=(r2(3)-r1(3))/ddist
 985 #ifdef DEBUG
 986       print *,"norm",er(1)**2+er(2)**2+er(3)**2
 987       print *,"r1",r1," r2",r2
 988       print *," er",er
 989       print *," dd",dd," dd3",dd3
 990 #endif
 991       do i=1,3
 992         do j=1,3
 993           a(i,j)=-3*er(i)*er(j)*dd3
 994 c          a(i,j)=-3*er(i)*er(j)
 995         enddo
 996       enddo
 997       do i=1,3
 998         a(i,i)=dd3+a(i,i)
 999 c        a(i,i)=1.0d0+a(i,i)
1000       enddo
1001 #ifdef DEBUG
1002       print *,"a before transformation",a
1003       print *,"emat1t",emat1t
1004       print *,"emat2",emat2
1005 #endif
1006       call matmat(emat1t,a,atemp)
1007 #ifdef DEBUG
1008       print *,"a after transformation",atemp
1009 #endif
1010       call matmat(atemp,emat2,a)
1011 #ifdef DEBUG
1012       print *,"a after transformation",a
1013       print *,"a and mu1"
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)
1018 #endif
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))
1021 #ifdef DEBUG
1022       print *,"eello_turn3_1",eello_turn3
1023 #endif
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
1040       do i=1,2
1041         do j=1,2
1042           e2(i,j)=e2gam(i,j)
1043         enddo
1044       enddo
1045 c      e1(1,1)=-e1(1,1)
1046 c      e1(1,2)=-e1(1,2)
1047       e2(1,1)=-e2(1,1)
1048       e2(1,2)=-e2(1,2)
1049       a(2,3)=-a(2,3)
1050       a(3,2)=-a(3,2)
1051       aux=0.0d0
1052       do i=1,2
1053         do j=1,2
1054           aux=aux+a(i+1,j+1)*(e1(i,1)*e2(1,j)+e1(i,2)*e2(2,j))
1055         enddo
1056       enddo
1057 #ifdef DEBUG
1058       print *,"eello_turn3_2",aux/2
1059 #endif
1060       eello_turn3=eello_turn3+aux/2
1061 c      stop
1062 c Derivatives in the elements of complete e1 and e2 matrices
1063       do i=1,2
1064         do j=1,2
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
1067         enddo
1068       enddo
1069 c Derivatives in the components of e1
1070       do i=1,2
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
1075       enddo
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
1081       do i=1,2
1082        eello_turn3_e2(1,1,i)=-(etemp2(1,i)*cosg+etemp2(2,i)*sing)
1083      &   *sint2sq
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)
1086      &   *sint2sq
1087        eello_turn3_e2(2,2,i)=eello_turn3_e2(1,2,i)*cost2
1088       enddo
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
1096       return
1097       end
1098 c-----------------------------------------------------------------------------
1099       subroutine edipole(mui,muj,emati,ematj,facpp,ene,enebi,enebj,
1100      &  idiri,idirj,lder)
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
1106       integer idiri,idirj
1107       logical lder
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
1120       else
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
1129       endif
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
1139       else
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
1148       endif
1149       return
1150       end
1151 c-----------------------------------------------------------------------------
1152       subroutine matvecsub(mat,vec)
1153       implicit none
1154       double precision mat(3,3),vec(3),auxvec(3)
1155       integer i,j
1156       do i=1,3
1157         auxvec(i)=0.0d0
1158         do j=1,3
1159           auxvec(i)=auxvec(i)+mat(i,j)*vec(j)
1160         enddo
1161       enddo
1162       vec=auxvec
1163       return
1164       end
1165 c-------------------------------------------------------------------------------
1166       subroutine rotmatsub(mat,gam,matgam)
1167       implicit none
1168       double precision mat(3,3),gam,matgam(3,3),cosg,sing
1169       integer i
1170       cosg=dcos(gam)
1171       sing=dsin(gam)
1172       do i=1,3
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
1176       enddo
1177       return
1178       end
1179 c-----------------------------------------------------------------------------
1180       subroutine coorsys(theta,phi,emat)
1181       implicit none
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
1188       cost=dcos(theta)
1189       sint=dsin(theta)
1190       cosphi=dcos(phi)
1191       sinphi=dsin(phi)
1192 #ifdef DEBUG
1193       print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1194 #endif
1195       emat=0.0d0
1196       emat(1,1)=sint*cosphi
1197       emat(2,1)=sint*sinphi
1198       emat(3,1)=cost
1199       emat(1,2)=-sinphi
1200       emat(2,2)=cosphi
1201       emat(3,2)=0.0d0
1202       emat(1,3)=-cost*cosphi
1203       emat(2,3)=-cost*sinphi
1204       emat(3,3)=sint
1205       return
1206       end
1207 c-----------------------------------------------------------------------------
1208       subroutine coorsys1(theta,phi,emat)
1209       implicit none
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
1216       cost=dcos(theta)
1217       sint=dsin(theta)
1218       cosphi=dcos(phi)
1219       sinphi=dsin(phi)
1220 #ifdef DEBUG
1221       print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1222 #endif
1223       emat=0.0d0
1224       emat(1,1)=sint*cosphi
1225       emat(2,1)=cost
1226       emat(3,1)=-sint*sinphi
1227       emat(1,2)=-cost*cosphi
1228       emat(2,2)=sint
1229       emat(3,2)=cost*sinphi
1230       emat(1,3)=sinphi
1231       emat(2,3)=0.0d0
1232       emat(3,3)=cosphi
1233       return
1234       end
1235 c-------------------------------------------------------------------------------
1236       subroutine matmat(a,b,c)
1237       implicit none
1238       double precision a(3,3),b(3,3),c(3,3)
1239       integer i,j,k
1240       do i=1,3
1241         do j=1,3
1242           c(i,j)=0.0d0
1243           do k=1,3
1244             c(i,j)=c(i,j)+a(i,k)*b(k,j)
1245           enddo
1246         enddo
1247       enddo
1248       return
1249       end