source/maxlik/src-Fmatch/fforcematch.F

   1       double precision function chisquare_force(n,x,g,ider)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include 'mpif.h'
   7       include "COMMON.MPI"
   8       integer ierror
   9 #endif
  10       include "COMMON.IOUNITS"
  11       include "COMMON.FMATCH"
  12       include "COMMON.CLASSES"
  13       include "COMMON.WEIGHTS"
  14       include "COMMON.VMCPAR"
  15       include "COMMON.PROTFILES"
  16       include "COMMON.MD"
  17       include "COMMON.DERIV"
  18       include "COMMON.INTERACT"
  19       include "COMMON.ENERGIES"
  20       include "COMMON.CHAIN"
  21       include "COMMON.NAMES"
  22       integer n
  23       double precision x(n),g(n),gg(max_paropt)
  24       integer ider
  25       logical lder
  26       integer i,ii,iii,ind,jj,j,k,kk,l,ll,iprot,nforce
  27       double precision aux,chisquare_(maxprot),
  28      & g_(max_paropt),fmean,fmean_,fsquare,fsquare_,fcov,fcov_
  29       integer ipass_conf,istart_conf,iend_conf
  30       double precision tcpu_ini,tcpu
  31       double precision energia (0:max_ene)
  32       double precision force_(3,maxres2,maxstr)
  33 #ifdef DEBUG
  34       write (iout,*) "Enter chisquare_forces, ider",ider
  35       write (iout,*) "x:",x
  36       write (iout,*) "ww",ww(:n_ene)
  37       call flush(iout)
  38 #endif
  39       calc_grad=.true.
  40 #ifdef MPI
  41 c      write (iout,*) "Before broadcast ider",ider,Master
  42 c      call flush(iout)
  43       call MPI_Bcast(ider,1,MPI_INTEGER,Master,Comm1,ierror)
  44 c      write (iout,*) "After broadcast ider ",ider
  45 c      call flush(iout)
  46 #endif
  47       if (ider.eq.0) return
  48       lder=ider.eq.2
  49
  50       if (lder) g=0.0d0
  51
  52       do iprot=1,nprot
  53
  54       chisquare(iprot)=0.0d0
  55
  56       gg = 0.0d0
  57
  58       fmean=0.0d0
  59       fsquare=0.0d0
  60       fcov=0.0d0
  61
  62       if (.not.fmatch(iprot)) cycle
  63
  64       call restore_molinfo(iprot)
  65
  66       if (.not. mod_other_params) then
  67         open (ientin,file=bforcefiles(iprot),status="old",
  68      &    form="unformatted",access="direct",recl=lenrec_forces(iprot))
  69 #ifdef DEBUG
  70         write (iout,*) "bforcefiles",bforcefiles(iprot)
  71 #endif
  72       else
  73         open (icbase,file=bprotfiles_MD(iprot),status="old",
  74      &    form="unformatted",access="direct",recl=lenrec_MD(iprot))
  75       endif
  76 c      write (iout,*) "chisquare: files opened"
  77 c      call flush(iout)
  78       ipass_conf=0
  79 #ifdef MPI
  80       do i=indstart_MD(me1,iprot),indend_MD(me1,iprot),maxstr_proc
  81         istart_conf=i
  82         iend_conf=min0(i+maxstr_proc-1,indend_MD(me1,iprot))
  83 #else
  84       do i=1,ntot_work_MD(iprot),maxstr_proc
  85         istart_conf=i
  86         iend_conf=min0(i+maxstr_proc-1,ntot_MD(iprot))
  87 #endif
  88 c        write (iout,*) "i",i," istart_conf",istart_conf," iend_conf",
  89 c     &    iend_conf
  90 c        call flush(iout)
  91         ii=0
  92 c
  93 c Read the chunk of conformations off a DA scratchfile.
  94 c
  95         if (.not. mod_other_params) then
  96 #ifdef DEBUG
  97           write (iout,*) "chisquare: calling daread_forces",istart_conf,
  98      &      iend_conf
  99 #endif
 100           call daread_forces(iprot,istart_conf,iend_conf)
 101           do iii=istart_conf,iend_conf
 102             ii=ii+1
 103             do k=1,nsite(iprot)
 104               do l=1,3
 105                 aux=0.0d0
 106                 do j=1,n_ene
 107                   aux=aux+forctb(j,l,k,ii,iprot)*ww(j)
 108                 enddo
 109                 force(l,k,iii,iprot)=aux
 110               enddo
 111             enddo
 112           enddo
 113         else
 114           call daread_MDcoords(iprot,istart_conf,iend_conf)
 115           do iii=istart_conf,iend_conf
 116             ii=ii+1
 117             call restore_MDcoords(iprot,ii)
 118             call int_from_cart1(.false.)
 119 #ifdef DEBUG
 120             write (iout,*) "Force: CG Coordinates conf",iii,ii
 121             do k=1,nres
 122               if (itype(k).eq.10 .or. itype(k).eq.ntyp1) then
 123                 write (iout,'(a4,i5,3f10.3)')
 124      &            restyp(itype(k)),k,(c(j,k),j=1,3)
 125                 else
 126                   write (iout,'(a4,i5,3f10.3,5x,3f10.3)')
 127      &            restyp(itype(k)),k,(c(j,k),j=1,3),(c(j,k+nres),j=1,3)
 128                endif
 129             enddo
 130             call int_from_cart1(.true.)
 131 #endif
 132             call zerograd
 133             call etotal(energia(0))
 134             call cartgrad
 135             do k=1,nres
 136               do l=1,3
 137                 force(l,k,iii,iprot)=-gcart(l,k)
 138                 do j=1,n_ene
 139                   forctb(j,l,k,ii,iprot)=-gcompon(j,l,k)
 140                 enddo
 141               enddo
 142             enddo
 143             kk=nres
 144             do k=nnt,nct
 145               if (itype(k).eq.10 .or. itype(k).eq.ntyp1) cycle
 146               kk=kk+1
 147               do l=1,3
 148                 force(l,kk,iii,iprot)=-gxcart(l,k)
 149                 do j=1,n_ene
 150                   forctb(j,l,kk,ii,iprot)=-gcomponx(j,l,k)
 151                 enddo
 152               enddo
 153             enddo
 154 #ifdef CHECKFORCE
 155             write (iout,*) "Checking forces from components"
 156             write (iout,'(2a5,3a10)')"coor","site","full",
 157      &          "from compon","diff"
 158             do k=1,nsite(iprot)
 159               do l=1,3
 160                 aux=0.0d0
 161                 do j=1,n_ene
 162                   aux=aux+forctb(j,l,k,ii,iprot)*ww(j)
 163                 enddo
 164                 write(iout,'(2i5,3f10.5)') l,k,force(l,k,iii,iprot),aux,
 165      &             (force(l,k,iii,iprot)-aux)/force(l,k,iii,iprot)
 166               enddo
 167             enddo
 168 #endif
 169           enddo
 170         endif
 171 c Compute chisquae components
 172         ii=0
 173         do iii=istart_conf,iend_conf
 174 #ifdef DEBUG
 175           write (iout,*) "Forces: conformation",iii,ii
 176           jj=nres
 177           do j=1,nres
 178             if (itype(j).eq.10 .or. itype(j).eq.ntyp1) then
 179               write (iout,'(a4,i5,3f10.5)')
 180      &          restyp(itype(j)),j,(force(k,j,iii,iprot),k=1,3)
 181             else
 182               jj=jj+1
 183               write (iout,'(a4,i5,3f10.5,5x,3f10.5)')
 184      &          restyp(itype(j)),j,(force(k,j,iii,iprot),k=1,3),
 185      &           (force(k,jj,iii,iprot),k=1,3)
 186             endif
 187           enddo
 188 #endif
 189           ii=ii+1
 190           do k=1,nsite(iprot)
 191             do l=1,3
 192               chisquare(iprot)=chisquare(iprot)+
 193      &          (CGforce_MD(l,k,iii,iprot)-force(l,k,iii,iprot))**2
 194               fmean=fmean+force(l,k,iii,iprot)
 195               fsquare=fsquare+force(l,k,iii,iprot)**2
 196               fcov=fcov+force(l,k,iii,iprot)*CGforce_MD(l,k,iii,iprot)
 197             enddo
 198           enddo
 199           if (lder) then
 200           ll=0
 201           do l=1,n_ene
 202             if (imask(l).gt.0) then
 203             ll=ll+1
 204             do j=1,nsite(iprot)
 205               do k=1,3
 206                 gg(ll)=gg(ll)+(CGforce_MD(k,j,iii,iprot)
 207      &            -force(k,j,iii,iprot))*forctb(l,k,j,ii,iprot)
 208 c                write (iout,*) "iii",iii," ii",ii,"l",l," ll",ll," j",j
 209 c                write (iout,*) " residual",
 210 c     &           CGforce_MD(k,j,iii,iprot)-force(k,j,iii,iprot),
 211 c     &           " forctab",forctb(l,k,j,ii,iprot)
 212 c                write (iout,*) "gg",gg(ll)
 213               enddo
 214             enddo
 215             endif
 216           enddo
 217           endif
 218         enddo
 219       enddo
 220
 221       if (lder) then
 222         ll=0
 223         do l=1,n_ene
 224          if (imask(l).gt.0) then
 225            ll=ll+1
 226            g(ll)=g(ll)+wforce(iprot)*gg(ll)/
 227      &       (3*ntot_work_MD(iprot)*nsite(iprot))
 228          endif
 229         enddo
 230 c        write (iout,*) "chisquare ll:",ll
 231       endif
 232 c      write (iout,*) "fmean",fmean," fsquare",fsquare
 233 #ifdef MPI
 234       fmean_=fmean
 235       call MPI_Reduce(fmean_,fmean,1,MPI_DOUBLE_PRECISION,
 236      &   MPI_SUM,Master,Comm1,ierror)
 237       fsquare_=fsquare
 238       call MPI_Reduce(fsquare_,fsquare,1,MPI_DOUBLE_PRECISION,
 239      &   MPI_SUM,Master,Comm1,ierror)
 240       fcov_=fcov
 241       call MPI_Reduce(fcov_,fcov,1,MPI_DOUBLE_PRECISION,
 242      &   MPI_SUM,Master,Comm1,ierror)
 243 #endif
 244       nforce = 3*ntot_work_MD(iprot)*nsite(iprot)
 245       fmean=fmean/nforce
 246 c      write (iout,*) "ntot_work",ntot_work_MD(iprot)," nforce",nforce
 247 c      write (iout,*) "fmean",fmean/nforce," fsquare",fsquare/nforce,
 248 c     & " variance",(fsquare-nforce*fmean)/nforce
 249       rcorr(iprot)=(fcov-nforce*meanforce(iprot)*fmean)
 250      &  /dsqrt(nforce*(fsquare-nforce*fmean*fmean)*varforce(iprot))
 251       enddo !v iprot
 252 #ifdef DEBUG
 253       do iprot=1,nprot
 254         write (iout,*) "Protein",iprot," forces"
 255         write (iout,"(4x,a,t55,a)") "Backbone","Sidechain"
 256 #ifdef MPI
 257         do iii=indstart_MD(me1,iprot),indend_MD(me1,iprot)
 258 #else
 259         do iii=1,ntot_work_MD(iprot)
 260 #endif
 261           write (iout,*) "Conformation",iii
 262           ii=nres
 263           do i=1,nres
 264             if (itype(i).eq.10 .or.itype(i).eq.ntyp1) then
 265               write (iout,'(a4,i5,3e15.5)') restyp(itype(i)),i,
 266      &          (CGforce_MD(j,i,iii,iprot),j=1,3)
 267               write (iout,'(9x,3e15.5)') (force(j,i,iii,iprot),j=1,3)
 268             else
 269               ii=ii+1
 270               write (iout,'(a4,i5,3e15.5,5x,3e15.5)')restyp(itype(i)),i,
 271      &          (CGforce_MD(j,i,iii,iprot),j=1,3),
 272      &          (CGforce_MD(j,ii,iii,iprot),j=1,3)
 273               write (iout,'(9x,3e15.5,5x,3e15.5)')
 274      &          (force(j,i,iii,iprot),j=1,3),
 275      &          (force(j,ii,iii,iprot),j=1,3)
 276             endif
 277           enddo
 278         enddo
 279       enddo
 280
 281       write(iout,*)"chisquare",chisquare(:nprot)
 282 #endif
 283 c Compute total chisquare and gradient
 284 #ifdef MPI
 285       chisquare_=chisquare
 286       call MPI_Reduce(chisquare_,chisquare,nprot,MPI_DOUBLE_PRECISION,
 287      &   MPI_SUM,Master,Comm1,ierror)
 288       if (lder) then
 289         gg=g
 290 c        write (iout,*) "g before reduce",g(:n)
 291         call MPI_Reduce(gg,g,n,MPI_DOUBLE_PRECISION,
 292      &   MPI_SUM,Master,Comm1,ierror)
 293 c        write (iout,*) "g after reduce",g(:n)
 294       endif
 295 #endif
 296 c      write (iout,*) "chisquare",chisquare(:nprot)
 297       aux=0.0d0
 298       do iprot=1,nprot
 299         if (.not.fmatch(iprot)) cycle
 300         chisquare(iprot)=chisquare(iprot)/
 301      &    (3*ntot_work_MD(iprot)*nsite(iprot))
 302         aux=aux+0.5d0*wforce(iprot)*chisquare(iprot)
 303       enddo
 304 #ifdef DEBUG
 305       write (iout,*) "Leave chisquare_forces"
 306       call flush(iout)
 307 #endif
 308 #ifdef MPI
 309       if (ider.eq.3) then
 310       do iprot=1,nprot
 311 #ifdef DEBUG
 312       write (iout,*) "fforce gather me1",me1," indstart",
 313      & indstart_MD(me1,iprot)," indend",indend_MD(me1,iprot)
 314 #endif
 315       do i=indstart_MD(me1,iprot),indend_MD(me1,iprot)
 316         do j=1,nsite(iprot)
 317           do k=1,3
 318             force_(k,j,i-indstart_MD(me1,iprot)+1)=force(k,j,i,iprot)
 319             force(k,j,i,iprot)=0.0d0
 320           enddo
 321         enddo
 322       enddo
 323 #ifdef DEBUG
 324       write (iout,*) "force_ scount_MD",scount_MD(me1,iprot)
 325       do i=1,scount_MD(me1,iprot)
 326         do j=1,nsite(iprot)
 327           write (iout,'(2i5,3f10.5)') i,j,(force_(k,j,i),k=1,3)
 328         enddo
 329       enddo
 330       write (iout,*) "scount_MD",scount_MD(:nprocs-1,iprot)
 331       write (iout,*) "idispl_MD",idispl_MD(:nprocs-1,iprot)
 332 #endif
 333       call MPI_Gatherv(force_(1,1,1),scount_MD(me1,iprot),MPI_FORCE,
 334      &   force(1,1,1,iprot),scount_MD(0,iprot),idispl_MD(0,iprot),
 335      &   MPI_FORCE,Master,Comm1, IERROR)
 336 #ifdef DEBUG
 337       write (iout,*) "Gather ierror",ierror
 338       write (iout,*) "force after gather"
 339       do i=1,ntot_MD(iprot)
 340         do j=1,nsite(iprot)
 341           write (iout,'(2i5,3f10.5)') i,j,(force(k,j,i,iprot),k=1,3)
 342         enddo
 343       enddo
 344 #endif
 345       enddo
 346       endif ! iderp.eq.3
 347 #endif
 348       chisquare_force=aux
 349       return
 350       end