source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/maxlikpdb.F

   1       double precision function maxlik_pdb(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),sumlik_(0:ntyp,0:ntyp),
  12      &  partf_,emin_,der_partf_(maxvar)
  13 #endif
  14       include "COMMON.IOUNITS"
  15       include "COMMON.PMFDERIV"
  16       include "COMMON.WEIGHTS"
  17       include "COMMON.TORSION"
  18       include "COMMON.FFIELD"
  19       include "COMMON.PDBSTAT"
  20       include "COMMON.GEO"
  21       double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
  22      &  b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0: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),enei,boltz,emin,
  26      & energia(36*36*72)
  27       double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
  28      & phi,cost1,cost2,sint1,sint2,
  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      & vol
  32       double precision conv /0.01745329251994329576d0/
  33       double precision r0pp(3)  /5.2739d0,5.4561d0,5.2261d0/
  34       double precision epspp(3) /1.6027d0,1.4879d0,0.0779d0/
  35       double precision dCACA(2) /3.784904d0,3.810180d0/
  36       double precision facp1(4),facpp(4),facturn3,ddist
  37       logical lder
  38       double precision x(n),g(n)
  39       integer ipoint,iipoint,ip,it1,it2,itp1,itp2,i,j,k,l,ii,jj!,irt1,irt2
  40       integer ider
  41       double precision sumlik(0:ntyp,0:ntyp),partf,der_partf(maxvar)
  42       integer itheta1,itheta2,igam
  43       vol = dlog((incr_theta*deg2rad)**2*incr_gam*deg2rad)
  44       delta = pi
  45 #ifdef MPI
  46 c nfun.eq.-1 : no broadcast
  47 c n.eq.-1 : stop function evaluation
  48 c nfun.ne.-1 .and. n.ne.-1 : master distribute variables to slaves, all do their share
  49 c of computations
  50       call MPI_Bcast(ider,1,MPI_LOGICAL,Master,MPI_COMM_WORLD,ierror)
  51       if (ider.le.0) return
  52 #endif
  53       lder=ider.gt.1
  54 c      write (iout,*) "MAXLIKPDB: ider",ider," lder",lder
  55 c      call flush(iout)
  56       IF (lder) THEN
  57 c Zero out gradient and hessian
  58       do i=1,n
  59         g(i)=0.0d0
  60       enddo
  61       ENDIF
  62
  63       sum=0.0d0
  64
  65       sumlik=0.0d0
  66       do i=1,npoint_pdb
  67 #ifdef DEBUG
  68       write (iout,*) "point",i," it1",it1," it2",it2
  69 #endif
  70       it1 = ipdbtyp(1,i)
  71       it2 = ipdbtyp(2,i)
  72       theta1=zz_pdb(1,i)
  73       theta2=zz_pdb(2,i)
  74       gam=(zz_pdb(3,i)+delta-pi)
  75       call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,jac,lder)
  76 #ifdef DEBUG
  77       write (2,*) "i",i," it1",it1," it2",it2," theta1",theta1,
  78      &  " theta2",theta2," gamma",gam," enei",enei
  79 #endif
  80       sumlik(it1,it2)=sumlik(it1,it2)+beta_pdb*enei*resol(i)
  81       IF (lder) THEN
  82 c Contributions to the gradient and the hessian
  83       do k=1,n
  84         g(k)=g(k)+beta_pdb*jac(k)*resol(i)
  85       enddo
  86       if (mask_beta_PDB.gt.0)
  87      &  g(mask_beta_PDB)=g(mask_beta_PDB)+enei*resol(i)
  88 #ifdef DEBUG
  89       print *,"i",i," jac",jac(:n)," resol",resol(i)," g",g(:n)
  90 #endif
  91       ENDIF ! lder
  92       enddo  ! i
  93 #ifdef MPI
  94       sumlik_=sumlik
  95       call MPI_Reduce(sumlik_,sumlik,(ntyp+1)*(ntortyp+1),
  96      &      MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,ierror)
  97       if (lder) then
  98         g_=g
  99         call MPI_Reduce(g_(1),g(1),n,MPI_DOUBLE_PRECISION,
 100      &   MPI_SUM,Master,Comm1,ierror)
 101 c        write (iout,*) "MAXLIKPDB: after reducing gradient"
 102 c        call flush(iout)
 103       endif
 104 #endif
 105 #ifdef DEBUG
 106       do it1=0,ntortyp-1
 107         do it2=0,ntortyp-1
 108           write (iout,*) "sumlik",it1,it2,sumlik(it1,it2)
 109         enddo
 110       enddo
 111 #endif
 112 c      sumlik=0
 113 c      if (lder) g=0
 114 #define CIPISZCZE
 115 #ifdef CIPISZCZE
 116 c Compute the partition functions for consecutive pairs of residue types
 117       do it1=0,ntortyp-1
 118         do it2=0,ntortyp-1
 119 #ifdef ENEOUT
 120           write (iout,'(a,2i2,2i10)')
 121      &      "BEGIN",it1,it2,istart_calka,iend_calka
 122 #endif
 123           ii=0
 124           do itheta1=70,160-incr_theta,incr_theta
 125             do itheta2=70,160-incr_theta,incr_theta
 126               do igam=-180,180-incr_gam,incr_gam
 127                 ii=ii+1
 128                 if (ii.ge.istart_calka .and. ii.le.iend_calka) then
 129                   theta1=(itheta1+0.5d0)*deg2rad
 130                   theta2=(itheta2+0.5d0)*deg2rad
 131                   gam=(igam+0.5d0)*deg2rad
 132                   call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
 133      &              jac,.false.)
 134                   energia(ii)=enei
 135 #ifdef ENEOUT
 136 c                  write (iout,'(3i5,1pe15.5)') itheta1,itheta2,igam,enei
 137 #endif
 138                 endif
 139               enddo
 140             enddo
 141           enddo
 142           emin=energia(istart_calka)
 143           do ii=istart_calka+1,iend_calka
 144             if (energia(ii).lt.emin) emin=energia(ii)
 145           enddo
 146 #ifdef MPI
 147           emin_=emin
 148           call MPI_AllReduce(emin_,emin,1,
 149      &         MPI_DOUBLE_PRECISION,MPI_MIN,Comm1,ierror)
 150 #endif
 151 c          emin=0.0d0
 152 c          write (2,*) it1,it2," emin",emin
 153           partf=0.0d0
 154           der_partf=0.0d0
 155           ii=0
 156           do itheta1=70,160-incr_theta,incr_theta
 157             do itheta2=70,160-incr_theta,incr_theta
 158               do igam=-180,180-incr_gam,incr_gam
 159                 ii=ii+1
 160                 if (ii.ge.istart_calka .and. ii.le.iend_calka) then
 161                   theta1=(itheta1+0.5d0)*deg2rad
 162                   theta2=(itheta2+0.5d0)*deg2rad
 163                   gam=(igam+0.5d0)*deg2rad
 164 c                    call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
 165 c     &              jac,lder)
 166                   enei=energia(ii)-emin
 167                   boltz = dsin(theta1)*dsin(theta2)*dexp(-beta_pdb*enei)
 168                   partf=partf+boltz
 169 c                  write (iout,*) "it1",it1," it2",it2,
 170 c     &              " theta1",theta1," theta2",theta2,
 171 c     &              " gamma",gam," boltz",boltz," enei",enei,
 172 c     &              " partf",partf
 173                   if (lder) then
 174                     call torval_e_der(n,it1,it2,theta1,theta2,gam,enei,
 175      &              jac,lder)
 176                     do i=1,n
 177                       der_partf(i)=der_partf(i)-boltz*jac(i)
 178                     enddo
 179                     if (mask_beta_PDB.gt.0)
 180      &              der_partf(mask_beta_PDB)=der_partf(mask_beta_PDB)-
 181      &               boltz*enei/beta_pdb
 182                   endif
 183                 endif
 184               enddo
 185             enddo
 186           enddo
 187 #ifdef ENEOUT
 188           write (iout,'(a,2i2)') "END",it1,it2
 189 #endif
 190 c          write (iout,*) "partf before REDUCE",it1,it2,partf
 191 #ifdef MPI
 192           partf_=partf
 193           call MPI_Reduce(partf_,partf,1,MPI_DOUBLE_PRECISION,
 194      &           MPI_SUM,Master,Comm1,ierror)
 195 c          write (iout,*) "Afte REDUCE: partf_",partf_," partf",partf
 196           if (lder) then
 197             der_partf_=der_partf
 198             call MPI_Reduce(der_partf_(1),der_partf(1),n,
 199      &         MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,ierror)
 200 c            write (iout,*) "MAXLIKPDB: after reducing gradient"
 201 c            call flush(iout)
 202           endif
 203           if (me.eq.Master) then
 204 #endif
 205           sumlik(it1,it2)=sumlik(it1,it2)+dlog(partf)-beta_pdb*emin+vol
 206 #ifdef DEBUG
 207           write(iout,*) "it1",it1," it2",it2," partf",partf," sumlik",
 208      &      sumlik(it1,it2)," emin",emin
 209 #endif
 210           if (lder) then
 211             do i=1,n
 212               g(i)=g(i)+beta_pdb*der_partf(i)/partf
 213             enddo
 214           endif
 215 #ifdef MPI
 216         endif
 217 #endif
 218         enddo
 219       enddo
 220 #ifdef ENEOUT
 221       stop
 222 #endif
 223 c      write (iout,*) "MAXLIKPDB: end of loop"
 224 c      call flush(iout)
 225 #endif
 226 c Compute the total contribution to the target function from a given core
 227       sum=0.0d0
 228       do it1=0,ntortyp-1
 229         do it2=0,ntortyp-1
 230           sum=sum+sumlik(it1,it2)
 231         enddo
 232       enddo
 233       maxlik_pdb = sum
 234       return
 235       end
 236 c----------------------------------------------------------------------
 237       subroutine torval_e_der(nvar,it1,it2,theta1,theta2,gam,enei,jac,
 238      &  lder)
 239       implicit none
 240       include "DIMENSIONS"
 241       include "DIMENSIONS.ZSCOPT"
 242       integer n,nmax
 243       parameter (nmax=max_paropt)
 244 #ifdef MPI
 245       include 'mpif.h'
 246       include 'COMMON.MPI'
 247       integer ierror
 248       double precision g_(max_paropt),sum_
 249 #endif
 250       include "COMMON.IOUNITS"
 251       include "COMMON.PMFDERIV"
 252       include "COMMON.WEIGHTS"
 253       include "COMMON.TORSION"
 254       include "COMMON.FFIELD"
 255       include "COMMON.GEO"
 256       integer nvar
 257       integer it1,it2,itp1,itp2,ip,i,j,k,l,ii,jj
 258       double precision b1m(3,2),b2m(3,2),dm(2,2),cm(2,2),
 259      &  b1m2(3,2),b2m1(3,2),b11,b12,b21,b22,c11,c12,d11,d12,c1(0:3),
 260      &  c2(0:3),em(2,2,2),em2(2,2,2),e0m(3),e0m2(3),jac(nmax)
 261       double precision v1_kcc_loc(3,3,2),v2_kcc_loc(3,3,2),
 262      & etoriv1(3,3,2),etoriv2(3,3,2),enei,etori,ebendi1,ebendi2
 263       double precision theta1,theta2,gam,thetap1,thetap2,gam1,gam2,
 264      & phi,cost1,cost2,sint1,sint2,
 265      & cosg,sing,cos2g,sin2g,sum,d3,delta,f,fb11,fb21,fb12,fb22,fc11,
 266      & fc12,fd11,fd12,sint1sint2,cosphi,sinphi,sumvalc,sumvals
 267       integer ind_ang /11/, ind_tor /13/
 268       logical lder
 269       double precision tschebyshev
 270 c following to test with selected types
 271       if (it1.le.1) then
 272         itp1=1
 273       else
 274         itp1=2
 275       endif
 276       if (iabs(it2).le.1) then
 277         itp2=1
 278       else
 279         itp2=2
 280       endif
 281
 282       ip=ip+1
 283
 284       delta=pi
 285
 286       do i=1,2
 287         do k=1,3
 288           b1m(k,i)=bnew1tor(k,i,it2)
 289           b1m2(k,i)=bnew2tor(k,i,it2)
 290           b2m1(k,i)=bnew1tor(k,i,it1)
 291           b2m(k,i)=bnew2tor(k,i,it1)
 292         enddo
 293       enddo
 294       do i=1,2
 295         do k=1,2
 296           cm(k,i)=ccnewtor(k,i,it2)
 297           dm(k,i)=ddnewtor(k,i,it1)
 298         enddo
 299       enddo
 300       do i=1,2
 301         do j=1,2
 302           do k=1,2
 303             em(k,j,i)=eenewtor(k,j,i,it1)
 304             em2(k,j,i)=eenewtor(k,j,i,it2)
 305           enddo
 306         enddo
 307       enddo
 308 c Invert coefficients for L-chirality
 309 #ifdef DEBUG
 310       write(iout,*) "it1",it1," it2",it2
 311       write(iout,*) "b1i (b1m)"
 312       do i=1,3
 313         write(iout,'(2f10.5,5x,2f10.5)') b1m(i,1),b1m(i,2)
 314       enddo
 315       write(iout,*) "b2j (b1m2)"
 316       do i=1,3
 317         write(iout,'(2f10.5,5x,2f10.5)') b1m2(i,1),b1m2(i,2)
 318       enddo
 319       write(iout,*) "b1j (b2m1)"
 320       do i=1,3
 321         write(iout,'(2f10.5,5x,2f10.5)') b2m1(i,1),b2m1(i,2)
 322       enddo
 323       write(iout,*) "b2j (b2m)"
 324       do i=1,3
 325         write(iout,'(2f10.5,5x,2f10.5)') b2m(i,1),b2m(i,2)
 326       enddo
 327       write (iout,*) "c"
 328       do i=1,2
 329         write (iout,'(2f10.5)') (cm(i,j),j=1,2)
 330       enddo
 331       write (iout,*) "d"
 332       do i=1,2
 333         write (iout,'(2f10.5)') (dm(i,j),j=1,2)
 334       enddo
 335 c      print *,"em1",em," em2",em2," e0m",e0m," e02m",e0m2
 336 #endif
 337       v1_kcc_loc=0.0d0
 338       v2_kcc_loc=0.0d0
 339       etoriv1=0.0d0
 340       etoriv2=0.0d0
 341       do k=1,3
 342         do l=1,3
 343           v1_kcc_loc(k,l,1)=b2m(l,1)*b1m(k,1)+b2m(l,2)*b1m(k,2)
 344           v2_kcc_loc(k,l,1)=b2m(l,2)*b1m(k,1)+b2m(l,1)*b1m(k,2)
 345         enddo
 346       enddo
 347       do k=1,2
 348         do l=1,2
 349           v1_kcc_loc(k,l,2)=-(cm(k,1)*dm(l,1)-cm(k,2)*dm(l,2))
 350           v2_kcc_loc(k,l,2)=-(cm(k,2)*dm(l,1)+cm(k,1)*dm(l,2))
 351         enddo
 352       enddo
 353       cost1=dcos(theta1)
 354       cost2=dcos(theta2)
 355       sint1=dsin(theta1)
 356       sint2=dsin(theta2)
 357       cosg=dcos(gam)
 358       sing=dsin(gam)
 359       cos2g=dcos(2*gam)
 360       sin2g=dsin(2*gam)
 361 c check with the formula from unres
 362       c1(0)=0.0d0
 363       c2(0)=0.0d0
 364       c1(1)=1.0d0
 365       c2(1)=1.0d0
 366       do j=2,3
 367         c1(j)=c1(j-1)*cost1
 368         c2(j)=c2(j-1)*cost2
 369       enddo
 370       etori=0.0d0
 371       sint1sint2=1.0d0
 372       do j=1,2
 373         cosphi=dcos(j*gam)
 374         sinphi=dsin(j*gam)
 375         sint1sint2=sint1sint2*sint1*sint2
 376         sumvalc=0.0d0
 377         sumvals=0.0d0
 378         do k=1,3
 379           do l=1,3
 380             sumvalc=sumvalc+v1_kcc_loc(l,k,j)*c1(k)*c2(l)
 381             sumvals=sumvals+v2_kcc_loc(l,k,j)*c1(k)*c2(l)
 382             etoriv1(l,k,j)=c1(k)*c2(l)*sint1sint2*cosphi
 383             etoriv2(l,k,j)=c1(k)*c2(l)*sint1sint2*sinphi
 384           enddo
 385         enddo
 386         etori=etori+(sumvalc*cosphi+sumvals*sinphi)*sint1sint2
 387       enddo
 388 c Bending contributions
 389       ebendi1=v1bend_chyb(0,it1)+
 390      &    tschebyshev(1,nbend_kcc_Tb(it1),v1bend_chyb(1,it1),cost1)
 391       ebendi2=v1bend_chyb(0,it2)+
 392      &    tschebyshev(1,nbend_kcc_Tb(it2),v1bend_chyb(1,it2),cost2)
 393 c      write (iout,*) "etori",etori," ebendi1",ebendi1," ebendi2",ebendi2
 394 c Compute derivatives of the energy in torsional parameters
 395       enei = wtor*etori+wang*(ebendi1+ebendi2)
 396 #ifdef ENEOUT
 397       write (iout,'(3f8.1,4(1pe15.5))') theta1*rad2deg,theta2*rad2deg,
 398      &    gam*rad2deg,etori,ebendi1,ebendi2,enei
 399 #endif
 400 c
 401       IF (lder) THEN
 402 c
 403       do k=1,nvar
 404         jac(k)=0.0d0
 405       enddo
 406 c Free term
 407 c Derivatives in b1 of the second residue
 408 c        jj=-1
 409       do j=1,3
 410 c          jj=jj+2
 411         jj=iabs(mask_bnew1tor(j,1,iabs(it2)))
 412         if (jj.gt.0) then
 413           do k=1,3
 414             jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,1)
 415      &        +etoriv2(j,k,1)*b2m(k,2)
 416           enddo
 417         endif
 418         jj=iabs(mask_bnew1tor(j,2,iabs(it2)))
 419         if (jj.gt.0) then
 420           do k=1,3
 421             if (it2.gt.0) then
 422               jac(jj)=jac(jj)+etoriv1(j,k,1)*b2m(k,2)
 423      &        +etoriv2(j,k,1)*b2m(k,1)
 424             else if (it2.lt.0) then
 425               jac(jj)=jac(jj)-etoriv1(j,k,1)*b2m(k,2)
 426      &        -etoriv2(j,k,1)*b2m(k,1)
 427             endif
 428           enddo
 429         endif
 430       enddo
 431 c Derivatives in b2 of the first residue
 432       do j=1,3
 433 c          jj=jj+2
 434         jj=iabs(mask_bnew2tor(j,1,it1))
 435         if (jj.gt.0) then
 436           do k=1,3
 437             jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,1)
 438      &      +etoriv2(k,j,1)*b1m(k,2)
 439           enddo
 440         endif
 441         jj=iabs(mask_bnew2tor(j,2,it1))
 442         if (jj.gt.0) then
 443           do k=1,3
 444             if (it1.gt.0) then
 445               jac(jj)=jac(jj)+etoriv1(k,j,1)*b1m(k,2)
 446      &        +etoriv2(k,j,1)*b1m(k,1)
 447             endif
 448           enddo
 449         endif
 450       enddo
 451 c Derivatives in c
 452       do j=1,2
 453 c        jj=jj+2
 454         jj=iabs(mask_ccnewtor(j,1,iabs(it2)))
 455         if (jj.gt.0) then
 456           do k=1,2
 457             jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,1)
 458      &       -etoriv2(j,k,2)*dm(k,2)
 459           enddo
 460         endif
 461         jj=iabs(mask_ccnewtor(j,2,iabs(it2)))
 462         if (jj.gt.0) then
 463           do k=1,2
 464             if (it2.gt.0) then
 465               jac(jj)=jac(jj)+etoriv1(j,k,2)*dm(k,2)
 466      &        -etoriv2(j,k,2)*dm(k,1)
 467             else if (it2.lt.0) then
 468               jac(jj)=jac(jj)-etoriv1(j,k,2)*dm(k,2)
 469      &        +etoriv2(j,k,2)*dm(k,1)
 470             endif
 471           enddo
 472         endif
 473       enddo
 474 c Derivatives in d
 475       do j=1,2
 476         jj=iabs(mask_ddnewtor(j,1,it1))
 477         if (jj.gt.0) then
 478           do k=1,2
 479             jac(jj)=jac(jj)-etoriv1(k,j,2)*cm(k,1)
 480      &       -etoriv2(k,j,2)*cm(k,2)
 481           enddo
 482         endif
 483         jj=iabs(mask_ddnewtor(j,2,it1))
 484         if (jj.gt.0) then
 485           do k=1,2
 486             if (it1.gt.0) then
 487               jac(jj)=jac(jj)+etoriv1(k,j,2)*cm(k,2)
 488      &        -etoriv2(k,j,2)*cm(k,1)
 489             endif
 490           enddo
 491         endif
 492       enddo
 493       do i=1,nvar
 494         jac(i)=wtor*jac(i)
 495       enddo
 496 c      jac=0.0d0
 497       do j=1,nbend_kcc_TB(it1)
 498         ii = mask_v1bend_chyb(j,it1)
 499         if (ii.gt.0) jac(ii) = wang*dcos(j*theta1)
 500       enddo
 501       do j=1,nbend_kcc_TB(it2)
 502         ii = mask_v1bend_chyb(j,it2)
 503         if (ii.gt.0) jac(ii) = jac(ii) + wang*dcos(j*theta2)
 504       enddo
 505       if (imask(ind_ang).gt.0) jac(imask(ind_ang))=ebendi1+ebendi2
 506       if (imask(ind_tor).gt.0) jac(imask(ind_tor))=etori
 507
 508       ENDIF
 509
 510       return
 511       end