source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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,2
  86       do it2=-2,2
  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.le.1) then
  95         itp1=1
  96       else
  97         itp1=2
  98       endif
  99       if (iabs(it2).le.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      else
 177 c        do i=1,2
 178 c          do k=1,3
 179 c            b1m(k,i)=bnew1(k,i,iabs(it2))
 180 c            b1m2(k,i)=bnew2(k,i,iabs(it2))
 181 c            b2m1(k,i)=bnew1(k,i,it1)
 182 c            b2m(k,i)=bnew2(k,i,it1)
 183 c          enddo
 184 c        enddo
 185 c        do i=1,2
 186 c          do k=1,2
 187 c            cm(k,i)=ccnew(k,i,iabs(it2))
 188 c            dm(k,i)=ddnew(k,i,it1)
 189 c          enddo
 190 c        enddo
 191 c        do i=1,2
 192 c          do j=1,2
 193 c            do k=1,3
 194 c              em(k,j,i)=eenew(k,j,i,it1)
 195 c              em2(k,j,i)=eenew(k,j,i,iabs(it2))
 196 c            enddo
 197 c          enddo
 198 c        enddo
 199 c        do k=1,3
 200 c          e0m(k)=e0new(k,it1)
 201 c          e0m2(k)=e0new(k,iabs(it2))
 202 c        enddo
 203 c       endif
 204 c      if (it2.gt.0) then
 205 c        b1m(1,1)=x(irt2+1)
 206 c        b1m(1,2)=x(irt2+2)
 207 c        b1m(2,1)=x(irt2+3)
 208 c        b1m(2,2)=x(irt2+4)
 209 c        b1m(3,1)=x(irt2+5)
 210 c        b1m(3,2)=x(irt2+6)
 211 c        b1m2(1,1)=x(irt2+7)
 212 c        b1m2(1,2)=x(irt2+8)
 213 c        b1m2(2,1)=x(irt2+9)
 214 c        b1m2(2,2)=x(irt2+10)
 215 c        b1m2(3,1)=x(irt2+11)
 216 c        b1m2(3,2)=x(irt2+12)
 217 c      else if (it2.eq.0) then
 218 c        b1m(1,1)=x(irt2+1)
 219 c        b1m(1,2)=0.0d0
 220 c        b1m(2,1)=x(irt2+3)
 221 c        b1m(2,2)=0.0d0
 222 c        b1m(3,1)=x(irt2+5)
 223 c        b1m(3,2)=0.0d0
 224 c        b1m2(1,1)=x(irt2+7)
 225 c        b1m2(1,2)=0.0d0
 226 c        b1m2(2,1)=x(irt2+9)
 227 c        b1m2(2,2)=0.0d0
 228 c        b1m2(3,1)=x(irt2+11)
 229 c        b1m2(3,2)=0.0d0
 230 c      else
 231 c        b1m(1,1)=x(irt2+1)
 232 c        b1m(1,2)=-x(irt2+2)
 233 c        b1m(2,1)=x(irt2+3)
 234 c        b1m(2,2)=-x(irt2+4)
 235 c        b1m(3,1)=x(irt2+5)
 236 c        b1m(3,2)=-x(irt2+6)
 237 c        b1m2(1,1)=x(irt2+7)
 238 c        b1m2(1,2)=-x(irt2+8)
 239 c        b1m2(2,1)=x(irt2+9)
 240 c        b1m2(2,2)=-x(irt2+10)
 241 c        b1m2(3,1)=x(irt2+11)
 242 c        b1m2(3,2)=-x(irt2+12)
 243 c      endif
 244 c      if (it1.eq.0) then
 245 c        b2m1(1,1)=x(irt1+1)
 246 c        b2m1(1,2)=0.0d0
 247 c        b2m1(2,1)=x(irt1+3)
 248 c        b2m1(2,2)=0.0d0
 249 c        b2m1(3,1)=x(irt1+5)
 250 c        b2m1(3,2)=0.0d0
 251 c        b2m(1,1)=x(irt1+7)
 252 c        b2m(1,2)=0.0d0
 253 c        b2m(2,1)=x(irt1+9)
 254 c        b2m(2,2)=0.0d0
 255 c        b2m(3,1)=x(irt1+11)
 256 c        b2m(3,2)=0.0d0
 257 c      else
 258 c        b2m1(1,1)=x(irt1+1)
 259 c        b2m1(1,2)=x(irt1+2)
 260 c        b2m1(2,1)=x(irt1+3)
 261 c        b2m1(2,2)=x(irt1+4)
 262 c        b2m1(3,1)=x(irt1+5)
 263 c        b2m1(3,2)=x(irt1+6)
 264 c        b2m(1,1)=x(irt1+7)
 265 c        b2m(1,2)=x(irt1+8)
 266 c        b2m(2,1)=x(irt1+9)
 267 c        b2m(2,2)=x(irt1+10)
 268 c        b2m(3,1)=x(irt1+11)
 269 c        b2m(3,2)=x(irt1+12)
 270 c      endif
 271 c      if (it2.gt.0) then
 272 c        cm(1,1)=x(irt2+13)
 273 c        cm(1,2)=x(irt2+14)
 274 c        cm(2,1)=x(irt2+15)
 275 c        cm(2,2)=x(irt2+16)
 276 c      else if (it2.eq.0) then
 277 c        cm(1,1)=x(irt2+13)
 278 c        cm(1,2)=0.0d0
 279 c        cm(2,1)=x(irt2+15)
 280 c        cm(2,2)=0.0d0
 281 c      else
 282 c        cm(1,1)=x(irt2+13)
 283 c        cm(1,2)=-x(irt2+14)
 284 c        cm(2,1)=x(irt2+15)
 285 c        cm(2,2)=-x(irt2+16)
 286 c      endif
 287 c      if (it1.eq.0) then
 288 c        dm(1,1)=x(irt1+17)
 289 c        dm(1,2)=0.0d0
 290 c        dm(2,1)=x(irt1+19)
 291 c        dm(2,2)=0.0d0
 292 c      else
 293 c        dm(1,1)=x(irt1+17)
 294 c        dm(1,2)=x(irt1+18)
 295 c        dm(2,1)=x(irt1+19)
 296 c        dm(2,2)=x(irt1+20)
 297 c      endif
 298 c E-coeffcients
 299 c      em2(1,1,1)=x(irt2+21)
 300 c      em2(2,1,1)=x(irt2+25)
 301 c      em2(1,2,2)=x(irt2+24)
 302 c      em2(2,2,2)=x(irt2+28)
 303 c      e0m2(1)=x(irt2+29)
 304 c      if (it2.gt.0) then
 305 c        em2(1,1,2)=x(irt2+22)
 306 c        em2(2,1,2)=x(irt2+26)
 307 c        em2(1,2,1)=x(irt2+23)
 308 c        em2(2,2,1)=x(irt2+27)
 309 c        em2(1,2,2)=x(irt2+24)
 310 c        em2(2,2,2)=x(irt2+28)
 311 c        e0m2(2)=x(irt2+30)
 312 c        e0m2(3)=x(irt2+31)
 313 c      else if (it2.eq.0) then
 314 c        em2(1,1,2)=0.0d0
 315 c        em2(2,1,2)=0.0d0
 316 c        em2(1,2,1)=0.0d0
 317 c        em2(2,2,1)=0.0d0
 318 c        e0m2(2)=0.0d0
 319 c        e0m2(3)=0.0d0
 320 c      else
 321 c        em2(1,1,2)=-x(irt2+22)
 322 c        em2(2,1,2)=-x(irt2+26)
 323 c        em2(1,2,1)=-x(irt2+23)
 324 c        em2(2,2,1)=-x(irt2+27)
 325 c        em2(1,2,2)=x(irt2+24)
 326 c        em2(2,2,2)=x(irt2+28)
 327 c        e0m2(2)=-x(irt2+30)
 328 c        e0m2(3)=-x(irt2+31)
 329 c      endif
 330 c      em(1,1,1)=x(irt1+21)
 331 c      em(2,1,1)=x(irt1+25)
 332 c      em(1,2,2)=x(irt1+24)
 333 c      em(2,2,2)=x(irt1+28)
 334 c      e0m(1)=x(irt1+29)
 335 c      if (it1.gt.0) then
 336 c        em(1,1,2)=x(irt1+22)
 337 c        em(2,1,2)=x(irt1+26)
 338 c        em(1,2,1)=x(irt1+23)
 339 c        em(2,2,1)=x(irt1+27)
 340 c        em(1,2,2)=x(irt1+24)
 341 c        em(2,2,2)=x(irt1+28)
 342 c        e0m(2)=x(irt1+30)
 343 c        e0m(3)=x(irt1+31)
 344 c      else if (it1.eq.0) then
 345 c        em(1,1,2)=0.0d0
 346 c        em(2,1,2)=0.0d0
 347 c        em(1,2,1)=0.0d0
 348 c        em(2,2,1)=0.0d0
 349 c        e0m(2)=0.0d0
 350 c        e0m(3)=0.0d0
 351 c      endif
 352 c Invert coefficients for L-chirality
 353 #ifdef DEBUG
 354       write(iout,*) "it1",it1," it2",it2
 355       write(iout,*) "b1i (b1m)"
 356       do i=1,3
 357         write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
 358       enddo
 359       write(iout,*) "b2j (b1m2)"
 360       do i=1,3
 361         write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
 362       enddo
 363       write(iout,*) "b1j (b2m1)"
 364       do i=1,3
 365         write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
 366       enddo
 367       write(iout,*) "b2j (b2m)"
 368       do i=1,3
 369         write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
 370       enddo
 371       write (iout,*) "c"
 372       do i=1,2
 373         write (iout,'(2f10.5)') (cm(i,j),j=1,2)
 374       enddo
 375       write (iout,*) "d"
 376       do i=1,2
 377         write (iout,'(2f10.5)') (dm(i,j),j=1,2)
 378       enddo
 379       print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
 380 #endif
 381       v1_kcc_loc=0.0d0
 382       v2_kcc_loc=0.0d0
 383       etoriv1=0.0d0
 384       etoriv2=0.0d0
 385 c      cm=0.0d0
 386 c      dm=0.0d0
 387       do k=1,3
 388         do l=1,3
 389           v1_kcc_loc(k,l,1)=b2m(l,1)*b1m(k,1)+b2m(l,2)*b1m(k,2)
 390           v2_kcc_loc(k,l,1)=b2m(l,2)*b1m(k,1)+b2m(l,1)*b1m(k,2)
 391         enddo
 392       enddo
 393       do k=1,2
 394         do l=1,2
 395           v1_kcc_loc(k,l,2)=-(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
 396           v2_kcc_loc(k,l,2)=-(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
 397         enddo
 398       enddo
 399 #define CHUJUN
 400 #ifdef CHUJUN
 401       if (torsion_pmf) then
 402
 403 c      write (iout,*) "it1",it1," it2",it2," etori",e0(it1,it2),
 404 c     &    " mask",mask_e0(it1,it2)
 405 c      write (iout,*)"mask_bnew1 1",it2,
 406 c     & (mask_bnew1tor(j,1,iabs(it2)),j=1,3)
 407 c      write (iout,*)"mask_bnew1 2",it2,
 408 c     & (mask_bnew1tor(j,2,iabs(it2)),j=1,3)
 409 c      write (iout,*)"mask_bnew2 1",it1,(mask_bnew2tor(j,1,it1),j=1,3)
 410 c      write (iout,*)"mask_bnew2 2",it1,(mask_bnew1tor(j,2,it1),j=1,3)
 411       iipoint=ipoint
 412       do ii=1,np(it1,it2)
 413         ipoint=ipoint+1
 414         i=ipoint
 415         theta1=zz(1,i)*conv
 416         theta2=zz(2,i)*conv
 417         gam=(zz(3,i)+delta-180)*conv
 418         cost1=dcos(theta1)
 419         cost2=dcos(theta2)
 420         sint1=dsin(theta1)
 421         sint2=dsin(theta2)
 422         cosg=dcos(gam)
 423         sing=dsin(gam)
 424         cos2g=dcos(2*gam)
 425         sin2g=dsin(2*gam)
 426 #ifdef OLDCALC
 427         b11=(b1m(1,1)+b1m(2,1)*cost2+b1m(3,1)*cost2*cost2)*sint2
 428         b12=(b1m(1,2)+b1m(2,2)*cost2+b1m(3,2)*cost2*cost2)*sint2
 429         b21=(b2m(1,1)+b2m(2,1)*cost1+b2m(3,1)*cost1*cost1)*sint1
 430         b22=(b2m(1,2)+b2m(2,2)*cost1+b2m(3,2)*cost1*cost1)*sint1
 431         d11=(dm(1,1)+dm(2,1)*cost1)*sint1*sint1
 432         d12=(dm(1,2)+dm(2,2)*cost1)*sint1*sint1
 433         c11=(cm(1,1)+cm(2,1)*cost2)*sint2*sint2
 434         c12=(cm(1,2)+cm(2,2)*cost2)*sint2*sint2
 435 #ifdef DEBUG
 436         write (iout,*) "b11",b11," b12",b12," b21",b21," b22",b22
 437         write (iout,*) "d11",d11," d12",d12," c11",c11," c12",c12
 438         write (iout,*) "cosg",cosg," cos2g",cos2g," sing",sing,
 439      &    " sin2g",sin2g
 440 #endif
 441         f=e0(it1,it2)+(b21*b11+b22*b12)*cosg+(b22*b11+b21*b12)*sing
 442      &    -(c11*d11-c12*d12)*cos2g-(c12*d11+c11*d12)*sin2g
 443         fb11=b21*cosg+b22*sing
 444         fb12=b22*cosg+b21*sing
 445         fb21=b11*cosg+b12*sing
 446         fb22=b12*cosg+b11*sing
 447         fc11=-d11*cos2g-d12*sin2g
 448         fc12= d12*cos2g-d11*sin2g
 449         fd11=-c11*cos2g-c12*sin2g
 450         fd12= c12*cos2g-c11*sin2g
 451 #endif
 452 c check with the formula from unres
 453         c1(0)=0.0d0
 454         c2(0)=0.0d0
 455         c1(1)=1.0d0
 456         c2(1)=1.0d0
 457         do j=2,3
 458           c1(j)=c1(j-1)*cost1
 459           c2(j)=c2(j-1)*cost2
 460         enddo
 461         etori=e0(it1,it2)
 462         sint1sint2=1.0d0
 463         do j=1,2
 464           cosphi=dcos(j*gam)
 465           sinphi=dsin(j*gam)
 466           sint1sint2=sint1sint2*sint1*sint2
 467           sumvalc=0.0d0
 468           sumvals=0.0d0
 469           do k=1,3
 470             do l=1,3
 471               sumvalc=sumvalc+v1_kcc_loc(l,k,j)*c1(k)*c2(l)
 472               sumvals=sumvals+v2_kcc_loc(l,k,j)*c1(k)*c2(l)
 473               etoriv1(l,k,j)=c1(k)*c2(l)*sint1sint2*cosphi
 474               etoriv2(l,k,j)=c1(k)*c2(l)*sint1sint2*sinphi
 475             enddo
 476           enddo
 477 c          if (j.eq.1) print *,"j",j," sumvlac",sumvalc*sint1sint2,
 478 c     &      b21*b11+b22*b12,
 479 c     &      "sumvals",sumvals*sint1sint2,b22*b11+b21*b12
 480 c          if (j.eq.2) print *,"j",j," sumvlac",sumvalc*sint1sint2,
 481 c     &      -(c11*d11-c12*d12),
 482 c     &      "sumvals",sumvals*sint1sint2,-(c12*d11+c11*d12)
 483           etori=etori+(sumvalc*cosphi+sumvals*sinphi)*sint1sint2
 484         enddo
 485 c      do k=1,3
 486 c        do l=1,3
 487 c          v1_kcc(k,l,1)=b1m(k,1)*b2m(l,1)+b1m(k,2)*b2m(l,2)
 488 c          v2_kcc(k,l,1)=b1m(k,1)*b2m(l,2)+b1m(k,2)*b2m(l,1)
 489 c        enddo
 490 c      enddo
 491 c      do k=1,2
 492 c        do l=1,2
 493 c          v1_kcc(k,l,2)=-0.25d0*(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
 494 c          v2_kcc(k,l,2)=-0.25d0*(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
 495 c        enddo
 496 c      enddo
 497 c        print *,"i",i," theta",zz(i,1),zz(i,2)," gamma",zz(i,3),
 498 c     &    " f",f," etori",etori
 499 #define CIPUN
 500 #ifdef CIPUN
 501         yc(1,i)=etori
 502 #ifdef DEBUG
 503         write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
 504      &    y(1,i),f,yc(1,i),y(1,i)-yc(1,i)
 505 #endif
 506         diff(1,i)=y(1,i)-yc(1,i)
 507         sum=sum+diff(1,i)*diff(1,i)*w(1,i)
 508 c
 509 c Compute derivatives of the energy in torsional parameters
 510 c
 511         IF (lder) THEN
 512 c
 513         do k=1,n
 514           jac(k)=0.0d0
 515         enddo
 516 c Free term
 517         if (mask_e0(it1,it2).gt.0) jac(mask_e0(it1,it2))=1.0d0
 518 c Derivatives in b1 of the second residue
 519 c        jj=-1
 520         do j=1,3
 521 c          jj=jj+2
 522           jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
 523           if (jj.gt.0) then
 524             do k=1,3
 525               jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
 526      &        +etoriv2(j,k,1)*b2m(k,2)
 527             enddo
 528           endif
 529           jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
 530           if (jj.gt.0) then
 531             do k=1,3
 532               if (it2.gt.0) then
 533                 jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,2)
 534      &          +etoriv2(j,k,1)*b2m(k,1)
 535               else if (it2.lt.0) then
 536                 jac(jj)=jac(jj)-etoriv1(j,k,1)*b2m(k,2)
 537      &          -etoriv2(j,k,1)*b2m(k,1)
 538               endif
 539             enddo
 540           endif
 541         enddo
 542 c Derivatives in b2 of the first residue
 543         do j=1,3
 544 c          jj=jj+2
 545           jj=iabs(mask_bnew2tor(j,1,it1))
 546           if (jj.gt.0) then
 547             do k=1,3
 548               jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
 549      &        +etoriv2(k,j,1)*b1m(k,2)
 550             enddo
 551           endif
 552           jj=iabs(mask_bnew2tor(j,2,it1))
 553           if (jj.gt.0) then
 554             do k=1,3
 555               if (it1.gt.0) then
 556                 jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
 557      &          +etoriv2(k,j,1)*b1m(k,1)
 558               endif
 559             enddo
 560           endif
 561         enddo
 562 c Derivatives in c
 563         do j=1,2
 564 c          jj=jj+2
 565           jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
 566           if (jj.gt.0) then
 567             do k=1,2
 568 c              jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,1)
 569 c     &         -0.25d0*etoriv2(j,k,2)*dm(k,2)
 570               jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,1)
 571      &         -etoriv2(j,k,2)*dm(k,2)
 572             enddo
 573           endif
 574           jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
 575           if (jj.gt.0) then
 576             do k=1,2
 577               if (it2.gt.0) then
 578 c                jac(jj+1)=jac(jj)+0.25d0*etoriv1(j,k,2)*dm(k,2)
 579 c     &          -0.25d0*etoriv2(j,k,2)*dm(k,1)
 580                 jac(jj)=jac(jj)+etoriv1(j,k,2)*dm(k,2)
 581      &          -etoriv2(j,k,2)*dm(k,1)
 582               else if (it2.lt.0) then
 583 c                jac(jj)=jac(jj)-0.25d0*etoriv1(j,k,2)*dm(k,2)
 584 c     &          +0.25d0*etoriv2(j,k,2)*dm(k,1)
 585                 jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,2)
 586      &          +etoriv2(j,k,2)*dm(k,1)
 587               endif
 588             enddo
 589           endif
 590         enddo
 591 c Derivatives in d
 592         do j=1,2
 593 c          jj=jj+2
 594           jj=iabs(mask_ddnewtor(j,1,it1))
 595           if (jj.gt.0) then
 596             do k=1,2
 597 c              jac(jj)=jac(jj)-0.25d0*etoriv1(k,j,2)*cm(k,1)
 598 c     &         -0.25d0*etoriv2(k,j,2)*cm(k,2)
 599               jac(jj)=jac(jj)-etoriv1(k,j,2)*cm(k,1)
 600      &         -etoriv2(k,j,2)*cm(k,2)
 601             enddo
 602           endif
 603           jj=iabs(mask_ddnewtor(j,2,it1))
 604           if (jj.gt.0) then
 605             do k=1,2
 606               if (it1.gt.0) then
 607 c                jac(jj)=jac(jj)+0.25d0*etoriv1(k,j,2)*cm(k,2)
 608 c     &          -0.25d0*etoriv2(k,j,2)*cm(k,1)
 609                 jac(jj)=jac(jj)+etoriv1(k,j,2)*cm(k,2)
 610      &          -etoriv2(k,j,2)*cm(k,1)
 611               endif
 612             enddo
 613           endif
 614         enddo
 615 c Contributions to the gradient and the hessian
 616         do k=1,n
 617           g(k)=g(k)+diff(1,i)*jac(k)*w(1,i)
 618 c          h(k,k)=h(k,k)+jac(k)**2*w(1,i)
 619 c          do l=1,k-1
 620 c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(1,i)
 621 c          enddo
 622         enddo
 623 #ifdef DEBUG
 624         print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
 625      &   " g",g(:n)
 626 #endif
 627         ENDIF ! lder
 628 #endif
 629       enddo  ! ii
 630
 631       endif ! torsion_pmf
 632
 633       do i=1,2
 634         do k=1,3
 635           b1m(k,i)=bnew1(k,i,it2)
 636           b1m2(k,i)=bnew2(k,i,it2)
 637           b2m1(k,i)=bnew1(k,i,it1)
 638           b2m(k,i)=bnew2(k,i,it1)
 639         enddo
 640       enddo
 641       do i=1,2
 642         do j=1,2
 643           do k=1,2
 644             em(k,j,i)=eenew(k,j,i,it1)
 645             em2(k,j,i)=eenew(k,j,i,it2)
 646           enddo
 647         enddo
 648       enddo
 649       do k=1,3
 650         e0m(k)=e0new(k,it1)
 651         e0m2(k)=e0new(k,it2)
 652       enddo
 653
 654       if (turn_pmf) then
 655
 656       ipoint=iipoint
 657       do ii=1,np(it1,it2)
 658         ipoint=ipoint+1
 659         i=ipoint
 660         theta1=zz(1,i)*conv
 661         theta2=zz(2,i)*conv
 662         gam=(zz(3,i)+delta-180)*conv
 663 #define CYCUN
 664 #ifdef CYCUN
 665 c eturn3 contributions
 666         call eturn3PMF(theta1,theta2,gam,b2m1,b1m2,em,em2,e0m,e0m2,
 667      &  facturn3,d1,d2,d3,eello_turn3,eello_turn3_b2,eello_turn3_b1,
 668      &  eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
 669         yc(2,i)=eello_turn3
 670 #ifdef DEBUG
 671         write (2,'(i7,3f7.1,4f10.3)') i,(zz(j,i),j=1,3),
 672      &    y(2,i),yc(2,i),y(2,i)-yc(2,i)
 673 #endif
 674         diff(2,i)=y(2,i)-yc(2,i)
 675         sum=sum+diff(2,i)*diff(2,i)*w(2,i)
 676
 677         IF (lder) THEN
 678
 679         jac=0.0d0
 680         if (mask_turn_PMF.gt.0) jac(mask_turn_PMF)=eello_turn3/wturn_PMF
 681 c Derivatives of in b1 of the first residue
 682 c        jj=-1
 683         do j=1,3
 684           jj=mask_bnew1(j,1,it1)
 685           if (jj.gt.0) then
 686 c            jj=jj+2
 687             jac(jj)=jac(jj)+eello_turn3_b2(j,1)
 688           endif
 689         enddo
 690         do j=1,3
 691           jj=mask_bnew1(j,2,it1)
 692           if (jj.gt.0) then
 693             if (it1.ne.0) then
 694               jac(jj)=jac(jj)+eello_turn3_b2(j,2)
 695             endif
 696           endif
 697         enddo
 698 c Derivatives of in b2 of the second residue
 699         do j=1,3
 700           jj=mask_bnew2(j,1,iabs(it2))
 701           if (jj.gt.0) then
 702             jac(jj)=jac(jj)+eello_turn3_b1(j,1)
 703           endif
 704         enddo
 705         do j=1,3
 706           jj=mask_bnew2(j,2,iabs(it2))
 707           if (jj.gt.0) then
 708             if (it2.gt.0) then
 709               jac(jj)=jac(jj)+eello_turn3_b1(j,2)
 710             else if (it2.lt.0) then
 711               jac(jj)=jac(jj)-eello_turn3_b1(j,2)
 712             endif
 713           endif
 714         enddo
 715 #define EE
 716 #ifdef EE
 717 c Derivatives in e of the first residue
 718 c        jj=21
 719         do j=1,2
 720           jj=mask_eenew(j,1,1,it1)
 721           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,1)
 722           jj=mask_eenew(j,2,2,it1)
 723           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,2)
 724           if (it1.ne.0) then
 725             jj=mask_eenew(j,1,2,it1)
 726             if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,1,2)
 727             jj=mask_eenew(j,2,1,it1)
 728             if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e1(j,2,1)
 729           endif
 730 c          jj=jj+4
 731         enddo
 732         jj=mask_e0new(1,it1)
 733         if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(1)
 734         if (it1.ne.0) then
 735           jj=mask_e0new(2,it1)
 736           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(2)
 737           jj=mask_e0new(3,it1)
 738           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m1(3)
 739         endif
 740 c Derivatives in e of the second residue
 741 c        jj=21
 742         do j=1,2
 743           jj=mask_eenew(j,1,1,iabs(it2))
 744           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,1)
 745           jj=mask_eenew(j,2,2,iabs(it2))
 746           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,2)
 747           if (it2.gt.0) then
 748             jj=mask_eenew(j,1,2,iabs(it2))
 749             if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,1,2)
 750             jj=mask_eenew(j,2,1,iabs(it2))
 751             if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e2(j,2,1)
 752           else if (it2.lt.0) then
 753             jj=mask_eenew(j,1,2,iabs(it2))
 754             if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,1,2)
 755             jj=mask_eenew(j,2,1,iabs(it2))
 756             if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e2(j,2,1)
 757           endif
 758 c          jj=jj+4
 759         enddo
 760         jj=mask_e0new(1,iabs(it2))
 761         if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(1)
 762         if (it2.gt.0) then
 763           jj=mask_e0new(2,iabs(it2))
 764           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(2)
 765           jj=mask_e0new(3,iabs(it2))
 766           if (jj.gt.0) jac(jj)=jac(jj)+eello_turn3_e0m2(3)
 767         else if (it2.lt.0) then
 768           jj=mask_e0new(2,iabs(it2))
 769           if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(2)
 770           jj=mask_e0new(3,iabs(it2))
 771           if (jj.gt.0) jac(jj)=jac(jj)-eello_turn3_e0m2(3)
 772         endif
 773 #endif
 774 c Contributions to the gradient and the hessian
 775         do k=1,n
 776           g(k)=g(k)+diff(2,i)*jac(k)*w(2,i)
 777 c          h(k,k)=h(k,k)+jac(k)**2*w(2,i)
 778 c          do l=1,k-1
 779 c            h(k,l)=h(k,l)+jac(k)*jac(l)*w(2,i)
 780 c          enddo
 781         enddo
 782 #ifdef DEBUG
 783         print *,"i",i," jac",jac(:n)," w, diff",w(1,i),diff(1,i),
 784      &   " g",g(:n)
 785 #endif
 786 c
 787         ENDIF ! lder
 788 #endif
 789       enddo
 790
 791       endif ! turn_pmf
 792 c      nfun=nfun+1
 793 c      rdif=sum
 794 c      return
 795 #endif
 796       if (eello_pmf) then
 797
 798       do ii=1,npel(it1,it2)
 799         eneelloc3b=0.0d0
 800         ipoint=ipoint+1
 801         theta1=zz(1,ipoint)*conv
 802         theta2=zz(2,ipoint)*conv
 803         ddist=zz(3,ipoint)**3
 804         thetap1=zz(4,ipoint)*conv
 805         thetap2=zz(5,ipoint)*conv
 806         phi=zz(6,ipoint)*conv
 807         gam1=zz(7,ipoint)*conv
 808         gam2=zz(8,ipoint)*conv
 809 #ifdef DEBUG
 810         print *,"ipoint",ipoint," theta1",theta1," theta2",theta2,
 811      &   " dd",dd," thetap1",thetap1," thetap2", thetap2," phi",phi,
 812      &   " gam1",gam1," gam2",gam2
 813 #endif
 814         do k=1,4
 815           facpp(k)=facp1(k)/ddist
 816         enddo
 817 #ifdef DEBUG
 818         print *,"facp1",facp1," facpp",facpp
 819 #endif
 820         call eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
 821      &    phi,b2m1,b2m,b1m,b1m2,facpp,eneelloc3,eneelloc3b,lder)
 822 #ifdef DEBUG
 823         print *,"eneelloc3",eneelloc3
 824 #endif
 825         do k=1,4
 826           yc(k,ipoint)=eneelloc3(k)
 827           diff(k,ipoint)=y(k,ipoint)-yc(k,ipoint)
 828           sum=sum+diff(k,ipoint)*diff(k,ipoint)*w(k,ipoint)
 829         enddo
 830 #ifdef DEBUG
 831         write (2,'(i7,7f7.1,4(2f10.3,5x))') i,(zz(j,ipoint),j=1,7),
 832      &    (y(k,ipoint),yc(k,ipoint),k=1,4)
 833 #endif
 834         IF (lder) THEN
 835         do l=1,4
 836 c          print *,"ipoint",ipoint," l",l
 837           do k=1,n
 838             jac(k)=0.0d0
 839           enddo
 840           if (mask_ello_PMF.gt.0)
 841      &    jac(mask_ello_PMF)=jac(mask_ello_PMF)+eneelloc3(l)/wello_PMF
 842 c#ifndef TORUN
 843 c          jac(16)=0.0d0
 844 c#endif
 845 c          jac(2)=eneelloc3(l)/wello
 846 c Derivatives in b1 of the first residue
 847 c          jj=-1
 848           do j=1,3
 849             jj=mask_bnew1(j,1,it1)
 850 c            jj=jj+2
 851             if (jj.gt.0) jac(jj)=eneelloc3b(j,1,1,1,l)
 852             jj=mask_bnew1(j,2,it1)
 853             if (jj.gt.0) then
 854               if (it1.gt.0) then
 855                 jac(jj)=eneelloc3b(j,2,1,1,l)
 856               else if (it1.lt.0) then
 857                 jac(jj)=-eneelloc3b(j,2,1,1,l)
 858               endif
 859             endif
 860           enddo
 861 c          print *,"bi1: jac",jac(:n)
 862 c Derivatives in b2 of the first residue
 863           do j=1,3
 864 c            jj=jj+2
 865             jj=mask_bnew2(j,1,it1)
 866             if (jj.gt.0) jac(jj)=eneelloc3b(j,1,2,1,l)
 867             jj=mask_bnew2(j,2,it1)
 868             if (jj.gt.0) then
 869               if (it1.gt.0) then
 870                 jac(jj)=eneelloc3b(j,2,2,1,l)
 871               else if (it1.lt.0) then
 872                 jac(jj)=-eneelloc3b(j,2,2,1,l)
 873               endif
 874             endif
 875           enddo
 876 c          print *,"bi2: jac",jac(:n)
 877 c Derivatives in b1 of the second residue
 878 c          jj=-1
 879           do j=1,3
 880             jj=mask_bnew1(j,1,iabs(it2))
 881             if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,1,2,l)
 882             jj=mask_bnew1(j,2,iabs(it2))
 883             if (jj.gt.0) then
 884               if (it2.gt.0) then
 885                 jac(jj)=jac(jj)+eneelloc3b(j,2,1,2,l)
 886               else if (it2.lt.0) then
 887                 jac(jj)=jac(jj)-eneelloc3b(j,2,1,2,l)
 888               endif
 889             endif
 890           enddo
 891 c          print *,"bj1: jac",jac(:n)
 892 c Derivatives in b2 of the second residue
 893           do j=1,3
 894             jj=mask_bnew2(j,1,iabs(it2))
 895             if (jj.gt.0) jac(jj)=jac(jj)+eneelloc3b(j,1,2,2,l)
 896             jj=mask_bnew2(j,2,iabs(it2))
 897             if (jj.gt.0) then
 898               if (it2.gt.0) then
 899                 jac(jj)=jac(jj)+eneelloc3b(j,2,2,2,l)
 900               else if (it2.lt.0) then
 901                 jac(jj)=jac(jj)-eneelloc3b(j,2,2,2,l)
 902               endif
 903             endif
 904           enddo
 905 c          print *,"bj2: jac",jac(:n)
 906           do k=1,n
 907             g(k)=g(k)+diff(l,ipoint)*jac(k)*w(l,ipoint)
 908 c            h(k,k)=h(k,k)+jac(k)**2*w(l,ipoint)
 909 c            do ll=1,k-1
 910 c              h(k,ll)=h(k,ll)+jac(k)*jac(ll)*w(l,ipoint)
 911 c            enddo
 912           enddo
 913 #ifdef DEBUG
 914           print *,"ipoint",ipoint," jac",jac(:n)," w, diff",w(l,ipoint),
 915      &     diff(1,ipoint)," g",g(:n)
 916 #endif
 917         enddo ! l
 918 c Contributions to the gradient and the hessian
 919
 920         ENDIF
 921       enddo ! ii
 922
 923       endif ! eello_PMF
 924
 925       enddo ! it2
 926       enddo ! it1
 927
 928 #ifdef MPI
 929 c      if (nfun.ne.-1) then
 930       sum_=sum
 931 c      write (2,*) "sum",sum
 932       call flush(2)
 933       call MPI_Reduce(sum_,sum,1,MPI_DOUBLE_PRECISION,
 934      &   MPI_SUM,Master,Comm1,ierror)
 935 c      endif
 936 c      write (iout,*) "after reducing function lder",lder
 937 c      call flush(iout)
 938       if (lder) then
 939         g_=g
 940         call MPI_Reduce(g_(1),g(1),n,MPI_DOUBLE_PRECISION,
 941      &   MPI_SUM,Master,Comm1,ierror)
 942 c      write (iout,*) "after reducing gradient"
 943 c      call flush(iout)
 944       endif
 945 #endif
 946
 947 c      if (nfun.ne.-1) nfun=nfun+1
 948
 949       rdif=0.5d0*sum
 950 c      write (iout,*) "Exit rdif"
 951 c      call flush(iout)
 952 c      stop
 953       return
 954       end
 955 c----------------------------------------------------------------------
 956       subroutine eello3(theta1,theta2,gam1,gam2,thetap1,thetap2,
 957      &    phi,bi1,bi2,bj1,bj2,facpp,eneelloc3,eneelloc3b,lder)
 958       implicit none
 959       include "DIMENSIONS"
 960       include "DIMENSIONS.ZSCOPT"
 961       include "COMMON.IOUNITS"
 962       double precision theta1,theta2,gam1,gam2,dd,thetap1,thetap2,phi
 963       double precision bi1(3,2),bi2(3,2),bj1(3,2),bj2(3,2),elpp(2,2),
 964      &    r0pp(2,2)
 965       double precision eneelloc3(4),eneelloc3b(3,2,2,2,4),emat1(3,3),
 966      &    emat2(3,3),emati1(3,3),emati2(3,3),ematj1(3,3),ematj2(3,3)
 967       double precision mui1(3),mui2(3),muj1(3),muj2(3)
 968       double precision cost1,cost2,sint1,sint2,costp1,sintp1
 969       double precision facpp(4)
 970       logical lder
 971       common /calcdip/ cost1,cost2,sint1,sint2
 972       double precision pi /3.14159265358979323844d0/
 973       integer j,k
 974       mui1=0.0d0
 975       mui2=0.0d0
 976       muj1=0.0d0
 977       muj2=0.0d0
 978       cost1=dcos(theta1)
 979       sint1=dsin(theta1)
 980       cost2=dcos(theta2)
 981       sint2=dsin(theta2)
 982 #ifdef DEBUG
 983       do j=1,3
 984       print '(4(a,2f10.5))',"bi1",(bi1(j,k),k=1,2)," bi2",
 985      &  (bi2(j,k),k=1,2),
 986      & " bj1",(bj1(j,k),k=1,2)," bj2",(bj2(j,k),k=1,2)
 987       enddo
 988       print *,"cost1",cost1," sint1",sint1," cost2",cost2," sint2",sint2
 989 #endif
 990       do j=1,2
 991         mui1(j+1)=(bi1(1,j)+bi1(2,j)*cost1+bi1(3,j)*cost1*cost1)*sint1
 992         muj1(j+1)=(bj1(1,j)+bj1(2,j)*cost2+bj1(3,j)*cost2*cost2)*sint2
 993         mui2(j+1)=(bi2(1,j)+bi2(2,j)*cost1+bi2(3,j)*cost1*cost1)*sint1
 994         muj2(j+1)=(bj2(1,j)+bj2(2,j)*cost2+bj2(3,j)*cost2*cost2)*sint2
 995       enddo
 996 c      mui2(2)=-mui2(2)
 997 c      muj2(2)=-muj2(2)
 998       mui1(2)=-mui1(2)
 999       muj1(2)=-muj1(2)
1000 #ifdef DEBUG
1001       print *,"Before rotation"
1002       print '(a,3f10.5)',"mui1",mui1
1003       print '(a,3f10.5)',"muj1",muj1
1004       print '(a,3f10.5)',"mui2",mui2
1005       print '(a,3f10.5)',"muj2",muj2
1006 #endif
1007       sintp1=dsin(thetap1)
1008       costp1=dcos(thetap1)
1009 c      print *,"costp1",costp1," sintp1",sintp1
1010       emat1(1,1)=sintp1
1011       emat1(1,2)=0.0d0
1012       emat1(1,3)=-costp1
1013       emat1(2,1)=0.0d0
1014       emat1(2,2)=1.0d0
1015       emat1(2,3)=0.0d0
1016       emat1(3,1)=costp1
1017       emat1(3,2)=0.0d0
1018       emat1(3,3)=sintp1
1019       call coorsys(thetap2,phi,emat2)
1020 #ifdef DEBUG
1021       print *,"emat1"
1022       print '(3f10.5)',((emat1(j,k),k=1,3),j=1,3)
1023       print *,"emat2"
1024       print '(3f10.5)',((emat2(j,k),k=1,3),j=1,3)
1025 #endif
1026       call rotmatsub(emat1,gam1,emati1)
1027       call rotmatsub(emat1,gam1,emati2)
1028 #ifdef DEBUG
1029       print *,"emati1"
1030       print '(3f10.5)',((emati1(k,j),k=1,3),j=1,3)
1031       print *,"emati2"
1032       print '(3f10.5)',((emati2(k,j),k=1,3),j=1,3)
1033 #endif
1034       call matvecsub(emati1,mui1)
1035       call matvecsub(emati2,mui2)
1036       call rotmatsub(emat2,gam2,ematj1)
1037       call rotmatsub(emat2,gam2,ematj2)
1038 #ifdef DEBUG
1039       print *,"ematj1"
1040       print '(3f10.5)',((ematj1(k,j),k=1,3),j=1,3)
1041       print *,"ematj2"
1042       print '(3f10.5)',((ematj2(k,j),k=1,3),j=1,3)
1043 #endif
1044       call matvecsub(ematj1,muj1)
1045       call matvecsub(ematj2,muj2)
1046 #ifdef DEBUG
1047       print *,"After rotation"
1048       print '(a,3f10.5)',"mui1",mui1
1049       print '(a,3f10.5)',"muj1",muj1
1050       print '(a,3f10.5)',"mui2",mui2
1051       print '(a,3f10.5)',"muj2",muj2
1052 #endif
1053       call edipole(mui1,muj1,emati1,ematj1,facpp(1),eneelloc3(1),
1054      &  eneelloc3b(1,1,1,1,1),eneelloc3b(1,1,1,2,1),-1,-1,lder)
1055       call edipole(mui2,muj1,emati2,ematj1,facpp(2),eneelloc3(2),
1056      &  eneelloc3b(1,1,2,1,2),eneelloc3b(1,1,1,2,2),1,-1,lder)
1057       call edipole(mui1,muj2,emati1,ematj2,facpp(3),eneelloc3(3),
1058      &  eneelloc3b(1,1,1,1,3),eneelloc3b(1,1,2,2,3),-1,1,lder)
1059       call edipole(mui2,muj2,emati2,ematj2,facpp(4),eneelloc3(4),
1060      &  eneelloc3b(1,1,2,1,4),eneelloc3b(1,1,2,2,4),1,1,lder)
1061       return
1062       end
1063 c-----------------------------------------------------------------------------
1064       subroutine eturn3PMF(theta1,theta2,gam,b1m,b2m,em1,em2,e0m1,e0m2,
1065      &  facturn,d1,d2,d3,eello_turn3,eello_turn3_b1,eello_turn3_b2,
1066      &  eello_turn3_e1,eello_turn3_e2,eello_turn3_e0m1,eello_turn3_e0m2)
1067       implicit none
1068       double precision theta1,theta2,gam,cost1,cost2,sint1,sint2,cosg,
1069      &  sing,sint1sq,sint2sq,aux
1070       double precision b1m(3,2),b2m(3,2),b1(2),b2(2),em1(2,2,2),
1071      &  em2(2,2,2),e0m1(3),e0m2(3),e1(2,2),e2(2,2),facturn,emat1(3,3),
1072      &  emat2(3,3),emat1t(3,3),mu1(2),mu2(2),e2gam(2,2),eello_turn3,
1073      &  eello_turn3_b2(3,2),eello_turn3_b1(3,2),eello_turn3_e1(2,2,2),
1074      &  eello_turn3_e2(2,2,2),eello_turn3_e0m1(3),eello_turn3_e0m2(3),
1075      &  r1(3),r2(3),er(3),a(3,3),atemp(3,3),etemp1(2,2),etemp2(2,2),
1076      &  d1,d2,d3,ddist,dd3
1077       double precision pi /3.14159265358979323844d0/
1078       integer i,j,k,l
1079 c diagnostics
1080 c      d1=3.8
1081 c      d2=3.8
1082 c      d3=3.8
1083 c end diagnostics
1084       cost1=dcos(theta1)
1085       cost2=dcos(theta2)
1086       sint1=dsin(theta1)
1087       sint1sq=sint1*sint1
1088       sint2=dsin(theta2)
1089       sint2sq=sint2*sint2
1090       cosg=dcos(gam)
1091       sing=dsin(gam)
1092       mu1(1)=(b1m(1,1)+b1m(2,1)*cost1+b1m(3,1)*cost1*cost1)*sint1
1093       mu1(2)=(b1m(1,2)+b1m(2,2)*cost1+b1m(3,2)*cost1*cost1)*sint1
1094       mu2(1)=(b2m(1,1)+b2m(2,1)*cost2+b2m(3,1)*cost2*cost2)*sint2
1095       mu2(2)=(b2m(1,2)+b2m(2,2)*cost2+b2m(3,2)*cost2*cost2)*sint2
1096       mu1(1)=-mu1(1)
1097       mu2(1)=-mu2(1)
1098       do k=1,2
1099         do l=1,2
1100           e1(k,l)=(em1(1,k,l)+em1(2,k,l)*cost1)*sint1sq
1101           e2(k,l)=(em2(1,k,l)+em2(2,k,l)*cost2)*sint2sq
1102         enddo
1103       enddo
1104       e1(1,1)=e1(1,1)+e0m1(1)*cost1
1105       e1(1,2)=e1(1,2)+e0m1(2)+e0m1(3)*cost1
1106       e1(2,1)=e1(2,1)+e0m1(2)*cost1+e0m1(3)
1107       e1(2,2)=e1(2,2)-e0m1(1)
1108       e2(1,1)=e2(1,1)+e0m2(1)*cost2
1109       e2(1,2)=e2(1,2)+e0m2(2)+e0m2(3)*cost2
1110       e2(2,1)=e2(2,1)+e0m2(2)*cost2+e0m2(3)
1111       e2(2,2)=e2(2,2)-e0m2(1)
1112 #ifdef DEBUG
1113       print *,"mu1",mu1," mu2",mu2," e1",e1," e2",e2
1114 #endif
1115 c Coordinate system of the first  dipole
1116       emat1(1,1)=sint1
1117       emat1(1,2)=cost1
1118       emat1(1,3)=0.0d0
1119       emat1(2,1)=-cost1
1120       emat1(2,2)=sint1
1121       emat1(2,3)=0.0d0
1122       emat1(3,1)=0.0d0
1123       emat1(3,2)=0.0d0
1124       emat1(3,3)=1.0d0
1125       do i=1,3
1126         do j=1,3
1127           emat1t(i,j)=emat1(j,i)
1128         enddo
1129       enddo
1130 #ifdef DEBUG
1131       print *,"theta1",theta1*180/pi," theta2",theta2*180/pi,
1132      &  " gam",gam*180/pi
1133 #endif
1134 c Calculate the reference system of the second unit
1135       call coorsys1(theta2,gam,emat2)
1136 c Dipole-dipole interaction matrix
1137 #ifdef DEBUG
1138       print *,"d1",d1," d2",d2," d3",d3
1139       print *,"emat1",emat1," emat2",emat2
1140 #endif
1141       r1(1)=-d1/2
1142       r1(2)=0.0d0
1143       r1(3)=0.0d0
1144       call matvecsub(emat1,r1)
1145       r2(1)=-d3/2
1146       r2(2)=0.0d0
1147       r2(3)=0.0d0
1148 #ifdef DEBUG
1149       print *," r2",r2
1150 #endif
1151       call matvecsub(emat2,r2)
1152 #ifdef DEBUG
1153       print *," r2",r2
1154 #endif
1155       r2(2)=r2(2)+d2
1156       ddist=dsqrt((r2(1)-r1(1))**2+(r2(2)-r1(2))**2+(r2(3)-r1(3))**2)
1157       dd3=facturn/ddist**3
1158       er(1)=(r2(1)-r1(1))/ddist
1159       er(2)=(r2(2)-r1(2))/ddist
1160       er(3)=(r2(3)-r1(3))/ddist
1161 #ifdef DEBUG
1162       print *,"norm",er(1)**2+er(2)**2+er(3)**2
1163       print *,"r1",r1," r2",r2
1164       print *," er",er
1165       print *," dd",dd," dd3",dd3
1166 #endif
1167       do i=1,3
1168         do j=1,3
1169           a(i,j)=-3*er(i)*er(j)*dd3
1170 c          a(i,j)=-3*er(i)*er(j)
1171         enddo
1172       enddo
1173       do i=1,3
1174         a(i,i)=dd3+a(i,i)
1175 c        a(i,i)=1.0d0+a(i,i)
1176       enddo
1177 #ifdef DEBUG
1178       print *,"a before transformation",a
1179       print *,"emat1t",emat1t
1180       print *,"emat2",emat2
1181 #endif
1182       call matmat(emat1t,a,atemp)
1183 #ifdef DEBUG
1184       print *,"a after transformation",atemp
1185 #endif
1186       call matmat(atemp,emat2,a)
1187 #ifdef DEBUG
1188       print *,"a after transformation",a
1189       print *,"a and mu1"
1190       print "(1h[,2f10.5,1h],5x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
1191      &  mu1(1),mu1(2),a(2,2),a(2,3),mu2(1)
1192       print "(27x,1h|,2f10.5,1h|,5x,1h|,f10.5,1h|)",
1193      &  a(3,2),a(3,3),mu2(2)
1194 #endif
1195       eello_turn3=(a(2,2)*mu1(1)*mu2(1)+a(2,3)*mu1(1)*mu2(2)
1196      &            +a(3,2)*mu1(2)*mu2(1)+a(3,3)*mu1(2)*mu2(2))
1197 #ifdef DEBUG
1198       print *,"eello_turn3_1",eello_turn3
1199 #endif
1200       eello_turn3_b1(1,1)=-(a(2,2)*mu2(1)+a(2,3)*mu2(2))*sint1
1201       eello_turn3_b1(2,1)=eello_turn3_b1(1,1)*cost1
1202       eello_turn3_b1(3,1)=eello_turn3_b1(2,1)*cost1
1203       eello_turn3_b1(1,2)=(a(3,2)*mu2(1)+a(3,3)*mu2(2))*sint1
1204       eello_turn3_b1(2,2)=eello_turn3_b1(1,2)*cost1
1205       eello_turn3_b1(3,2)=eello_turn3_b1(2,2)*cost1
1206       eello_turn3_b2(1,1)=-(a(2,2)*mu1(1)+a(3,2)*mu1(2))*sint2
1207       eello_turn3_b2(2,1)=eello_turn3_b2(1,1)*cost2
1208       eello_turn3_b2(3,1)=eello_turn3_b2(2,1)*cost2
1209       eello_turn3_b2(1,2)=(a(2,3)*mu1(1)+a(3,3)*mu1(2))*sint2
1210       eello_turn3_b2(2,2)=eello_turn3_b2(1,2)*cost2
1211       eello_turn3_b2(3,2)=eello_turn3_b2(2,2)*cost2
1212       e2gam(1,1)=e2(1,1)*cosg+e2(2,1)*sing
1213       e2gam(1,2)=e2(1,2)*cosg+e2(2,2)*sing
1214       e2gam(2,1)=-e2(1,1)*sing+e2(2,1)*cosg
1215       e2gam(2,2)=-e2(1,2)*sing+e2(2,2)*cosg
1216       do i=1,2
1217         do j=1,2
1218           e2(i,j)=e2gam(i,j)
1219         enddo
1220       enddo
1221 c      e1(1,1)=-e1(1,1)
1222 c      e1(1,2)=-e1(1,2)
1223       e2(1,1)=-e2(1,1)
1224       e2(1,2)=-e2(1,2)
1225       a(2,3)=-a(2,3)
1226       a(3,2)=-a(3,2)
1227       aux=0.0d0
1228       do i=1,2
1229         do j=1,2
1230           aux=aux+a(i+1,j+1)*(e1(i,1)*e2(1,j)+e1(i,2)*e2(2,j))
1231         enddo
1232       enddo
1233 #ifdef DEBUG
1234       print *,"eello_turn3_2",aux/2
1235 #endif
1236       eello_turn3=eello_turn3+aux/2
1237 c      stop
1238 c Derivatives in the elements of complete e1 and e2 matrices
1239       do i=1,2
1240         do j=1,2
1241           etemp1(i,j)=(a(i+1,2)*e2(j,1)+a(i+1,3)*e2(j,2))/2
1242           etemp2(i,j)=(a(2,j+1)*e1(1,i)+a(3,j+1)*e1(2,i))/2
1243         enddo
1244       enddo
1245 c Derivatives in the components of e1
1246       do i=1,2
1247         eello_turn3_e1(1,1,i)=etemp1(1,i)*sint1sq
1248         eello_turn3_e1(2,1,i)=eello_turn3_e1(1,1,i)*cost1
1249         eello_turn3_e1(1,2,i)= etemp1(2,i)*sint1sq
1250         eello_turn3_e1(2,2,i)=eello_turn3_e1(1,2,i)*cost1
1251       enddo
1252 c Derivatives in the components of e01m
1253       eello_turn3_e0m1(1)=etemp1(1,1)*cost1-etemp1(2,2)
1254       eello_turn3_e0m1(2)=etemp1(1,2)+cost1*etemp1(2,1)
1255       eello_turn3_e0m1(3)=etemp1(1,2)*cost1+etemp1(2,1)
1256 c Derivatives in the elements of e2
1257       do i=1,2
1258        eello_turn3_e2(1,1,i)=-(etemp2(1,i)*cosg+etemp2(2,i)*sing)
1259      &   *sint2sq
1260        eello_turn3_e2(2,1,i)=eello_turn3_e2(1,1,i)*cost2
1261        eello_turn3_e2(1,2,i)=(-etemp2(1,i)*sing+etemp2(2,i)*cosg)
1262      &   *sint2sq
1263        eello_turn3_e2(2,2,i)=eello_turn3_e2(1,2,i)*cost2
1264       enddo
1265 c Derivatives in the elements of e02m
1266       eello_turn3_e0m2(1)=-(etemp2(1,1)*cosg+etemp2(2,1)*sing)*cost2
1267      &                 +etemp2(1,2)*sing-etemp2(2,2)*cosg
1268       eello_turn3_e0m2(2)=-etemp2(1,2)*cosg-etemp2(2,2)*sing+
1269      &  (-etemp2(1,1)*sing+etemp2(2,1)*cosg)*cost2
1270       eello_turn3_e0m2(3)=(-etemp2(1,2)*cosg-etemp2(2,2)*sing)*cost2
1271      &  -etemp2(1,1)*sing+etemp2(2,1)*cosg
1272       return
1273       end
1274 c-----------------------------------------------------------------------------
1275       subroutine edipole(mui,muj,emati,ematj,facpp,ene,enebi,enebj,
1276      &  idiri,idirj,lder)
1277       include "DIMENSIONS"
1278       include "DIMENSIONS.ZSCOPT"
1279       include "COMMON.IOUNITS"
1280       double precision mui(3),muj(3),emati(3,3),ematj(3,3),
1281      &  facpp,ene,enebi(3,2),enebj(3,2),muimuj,cost1,sint1,cost2,sint2
1282       integer idiri,idirj
1283       logical lder
1284       common /calcdip/ cost1,cost2,sint1,sint2
1285       muimuj=mui(1)*muj(1)+mui(2)*muj(2)+mui(3)*muj(3)
1286       ene=facpp*(muimuj-3*mui(3)*muj(3))
1287       if (idiri.gt.0) then
1288         enebi(1,1)=(emati(1,2)*muj(1)+emati(2,2)*muj(2)
1289      &   +emati(3,2)*muj(3)-3*muj(3)*emati(3,2))*sint1*facpp
1290         enebi(2,1)=enebi(1,1)*cost1
1291         enebi(3,1)=enebi(2,1)*cost1
1292         enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
1293      &   +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
1294         enebi(2,2)=enebi(1,2)*cost1
1295         enebi(3,2)=enebi(2,2)*cost1
1296       else
1297         enebi(1,1)=(-emati(1,2)*muj(1)-emati(2,2)*muj(2)
1298      &   -emati(3,2)*muj(3)+3*muj(3)*emati(3,2))*sint1*facpp
1299         enebi(2,1)=enebi(1,1)*cost1
1300         enebi(3,1)=enebi(2,1)*cost1
1301         enebi(1,2)=(emati(1,3)*muj(1)+emati(2,3)*muj(2)
1302      &   +emati(3,3)*muj(3)-3*muj(3)*emati(3,3))*sint1*facpp
1303         enebi(2,2)=enebi(1,2)*cost1
1304         enebi(3,2)=enebi(2,2)*cost1
1305       endif
1306       if (idirj.gt.0) then
1307         enebj(1,1)=(ematj(1,2)*mui(1)+ematj(2,2)*mui(2)
1308      &   +ematj(3,2)*mui(3)-3*mui(3)*ematj(3,2))*sint2*facpp
1309         enebj(2,1)=enebj(1,1)*cost2
1310         enebj(3,1)=enebj(2,1)*cost2
1311         enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
1312      &   +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
1313         enebj(2,2)=enebj(1,2)*cost2
1314         enebj(3,2)=enebj(2,2)*cost2
1315       else
1316         enebj(1,1)=(-ematj(1,2)*mui(1)-ematj(2,2)*mui(2)
1317      &   -ematj(3,2)*mui(3)+3*mui(3)*ematj(3,2))*sint2*facpp
1318         enebj(2,1)=enebj(1,1)*cost2
1319         enebj(3,1)=enebj(2,1)*cost2
1320         enebj(1,2)=(ematj(1,3)*mui(1)+ematj(2,3)*mui(2)
1321      &   +ematj(3,3)*mui(3)-3*mui(3)*ematj(3,3))*sint2*facpp
1322         enebj(2,2)=enebj(1,2)*cost2
1323         enebj(3,2)=enebj(2,2)*cost2
1324       endif
1325       return
1326       end
1327 c-----------------------------------------------------------------------------
1328       subroutine matvecsub(mat,vec)
1329       implicit none
1330       double precision mat(3,3),vec(3),auxvec(3)
1331       integer i,j
1332       do i=1,3
1333         auxvec(i)=0.0d0
1334         do j=1,3
1335           auxvec(i)=auxvec(i)+mat(i,j)*vec(j)
1336         enddo
1337       enddo
1338       vec=auxvec
1339       return
1340       end
1341 c-------------------------------------------------------------------------------
1342       subroutine rotmatsub(mat,gam,matgam)
1343       implicit none
1344       double precision mat(3,3),gam,matgam(3,3),cosg,sing
1345       integer i
1346       cosg=dcos(gam)
1347       sing=dsin(gam)
1348       do i=1,3
1349         matgam(i,1)=mat(i,1)
1350         matgam(i,2)= mat(i,2)*cosg+mat(i,3)*sing
1351         matgam(i,3)=-mat(i,2)*sing+mat(i,3)*cosg
1352       enddo
1353       return
1354       end
1355 c-----------------------------------------------------------------------------
1356       subroutine coorsys(theta,phi,emat)
1357       implicit none
1358       include "DIMENSIONS"
1359       include "DIMENSIONS.ZSCOPT"
1360       include "COMMON.IOUNITS"
1361       double precision theta,phi,emat(3,3)
1362       double precision cost,sint,cosphi,sinphi
1363 c      print *,"theta",theta," phi",phi
1364       cost=dcos(theta)
1365       sint=dsin(theta)
1366       cosphi=dcos(phi)
1367       sinphi=dsin(phi)
1368 #ifdef DEBUG
1369       print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1370 #endif
1371       emat=0.0d0
1372       emat(1,1)=sint*cosphi
1373       emat(2,1)=sint*sinphi
1374       emat(3,1)=cost
1375       emat(1,2)=-sinphi
1376       emat(2,2)=cosphi
1377       emat(3,2)=0.0d0
1378       emat(1,3)=-cost*cosphi
1379       emat(2,3)=-cost*sinphi
1380       emat(3,3)=sint
1381       return
1382       end
1383 c-----------------------------------------------------------------------------
1384       subroutine coorsys1(theta,phi,emat)
1385       implicit none
1386       include "DIMENSIONS"
1387       include "DIMENSIONS.ZSCOPT"
1388       include "COMMON.IOUNITS"
1389       double precision theta,phi,emat(3,3)
1390       double precision cost,sint,cosphi,sinphi
1391 c      print *,"theta",theta," phi",phi
1392       cost=dcos(theta)
1393       sint=dsin(theta)
1394       cosphi=dcos(phi)
1395       sinphi=dsin(phi)
1396 #ifdef DEBUG
1397       print *,"cost",cost," sint",sint," cosphi",cosphi," sinphi",sinphi
1398 #endif
1399       emat=0.0d0
1400       emat(1,1)=sint*cosphi
1401       emat(2,1)=cost
1402       emat(3,1)=-sint*sinphi
1403       emat(1,2)=-cost*cosphi
1404       emat(2,2)=sint
1405       emat(3,2)=cost*sinphi
1406       emat(1,3)=sinphi
1407       emat(2,3)=0.0d0
1408       emat(3,3)=cosphi
1409       return
1410       end
1411 c-------------------------------------------------------------------------------
1412       subroutine matmat(a,b,c)
1413       implicit none
1414       double precision a(3,3),b(3,3),c(3,3)
1415       integer i,j,k
1416       do i=1,3
1417         do j=1,3
1418           c(i,j)=0.0d0
1419           do k=1,3
1420             c(i,j)=c(i,j)+a(i,k)*b(k,j)
1421           enddo
1422         enddo
1423       enddo
1424       return
1425       end