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