source/wham/src-M/enecalc1.F

   1       subroutine enecalc(islice,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "DIMENSIONS.FREE"
   6 #ifdef MPI
   7       include "mpif.h"
   8       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
   9       include "COMMON.MPI"
  10 #endif
  11       include "COMMON.CHAIN"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.PROTFILES"
  14       include "COMMON.NAMES"
  15       include "COMMON.VAR"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.GEO"
  18       include "COMMON.FFIELD"
  19       include "COMMON.ENEPS"
  20       include "COMMON.LOCAL"
  21       include "COMMON.WEIGHTS"
  22       include "COMMON.INTERACT"
  23       include "COMMON.FREE"
  24       include "COMMON.ENERGIES"
  25       include "COMMON.CONTROL"
  26       include "COMMON.TORCNSTR"
  27       character*64 nazwa
  28       character*80 bxname
  29       character*3 liczba
  30       double precision qwolynes
  31       external qwolynes
  32       integer errmsg_count,maxerrmsg_count /100/
  33       double precision rmsnat,gyrate
  34       external rmsnat,gyrate
  35       double precision tole /1.0d-1/
  36       integer i,itj,ii,iii,j,k,l,licz
  37       integer ir,ib,ipar,iparm
  38       integer iscor,islice
  39       real*4 csingle(3,maxres2)
  40       double precision energ
  41       double precision temp
  42       integer ilen,iroof
  43       external ilen,iroof
  44       double precision energia(0:max_ene),rmsdev,efree,eini
  45       double precision fT(6),quot,quotl,kfacl,kfac /2.4d0/,T0 /3.0d2/
  46       double precision tt
  47       integer snk_p(MaxR,MaxT_h,Max_parm)
  48       logical lerr
  49       character*64 bprotfile_temp
  50       call opentmp(islice,ientout,bprotfile_temp)
  51       iii=0
  52       ii=0
  53       errmsg_count=0
  54       write (iout,*) "enecalc: nparmset ",nparmset
  55 #ifdef MPI
  56       do iparm=1,nParmSet
  57         do ib=1,nT_h(iparm)
  58           do i=1,nR(ib,iparm)
  59             snk_p(i,ib,iparm)=0
  60           enddo
  61         enddo
  62       enddo
  63       do i=indstart(me1),indend(me1)
  64 #else
  65       do iparm=1,nParmSet
  66         do ib=1,nT_h(iparm)
  67           do i=1,nR(ib,iparm)
  68             snk(i,ib,iparm)=0
  69           enddo
  70         enddo
  71       enddo
  72       do i=1,ntot
  73 #endif
  74         read(ientout,rec=i,err=101)
  75      &    ((csingle(l,k),l=1,3),k=1,nres),
  76      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
  77      &    nss,(ihpb(k),jhpb(k),k=1,nss),
  78      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
  79          if (indpdb.gt.0) then
  80            do k=1,nres
  81              do l=1,3
  82                c(l,k)=csingle(l,k)
  83              enddo
  84            enddo
  85            do k=nnt,nct
  86              do l=1,3
  87                c(l,k+nres)=csingle(l,k+nres)
  88              enddo
  89            enddo
  90            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
  91            q(nQ+1,iii+1)=rmsnat(iii+1)
  92          endif
  93          q(nQ+2,iii+1)=gyrate(iii+1)
  94 c        fT=T0*beta_h(ib,ipar)*1.987D-3
  95 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
  96         if (rescale_mode.eq.1) then
  97           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
  98 #if defined(FUNCTH)
  99           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 100           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 101 #elif defined(FUNCT)
 102           ft(6)=quot
 103 #else
 104           ft(6)=1.0d0
 105 #endif
 106           quotl=1.0d0
 107           kfacl=1.0d0
 108           do l=1,5
 109             quotl=quotl*quot
 110             kfacl=kfacl*kfac
 111             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 112           enddo
 113         else if (rescale_mode.eq.2) then
 114           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 115 #if defined(FUNCTH)
 116           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 117           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 118 #elif defined(FUNCT)
 119           ft(6)=quot
 120 #else
 121           ft(6)=1.0d0
 122 #endif
 123           quotl=1.0d0
 124           do l=1,5
 125             quotl=quotl*quot
 126             fT(l)=1.12692801104297249644d0/
 127      &         dlog(dexp(quotl)+dexp(-quotl))
 128           enddo
 129         else if (rescale_mode.eq.0) then
 130           do l=1,5
 131             fT(l)=1.0d0
 132           enddo
 133         else
 134           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
 135      &     rescale_mode
 136           call flush(iout)
 137           return1
 138         endif
 139
 140 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 141 c     &   " kfac",kfac,"quot",quot," fT",fT
 142         do j=1,2*nres
 143           do k=1,3
 144             c(k,j)=csingle(k,j)
 145           enddo
 146         enddo
 147         call int_from_cart1(.false.)
 148         ii=ii+1
 149         do iparm=1,nparmset
 150
 151         call restore_parm(iparm)
 152 #ifdef DEBUG
 153             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 154      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 155      &      wtor_d,wsccor,wbond
 156 #endif
 157         call etotal(energia(0),fT)
 158 #ifdef DEBUG
 159         write (iout,*) "Conformation",i
 160           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 161           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 162         call enerprint(energia(0),fT)
 163         write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 164         write (iout,*) "ftors",ftors
 165         call briefout(i,energia(0))
 166         temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
 167         write (iout,*) "temp", temp
 168         call pdbout(i,temp,energia(0),energia(0),0.0d0,0.0d0)
 169 #endif
 170         if (energia(0).ge.1.0d20) then
 171           write (iout,*) "NaNs detected in some of the energy",
 172      &     " components for conformation",ii+1
 173           write (iout,*) "The Cartesian geometry is:"
 174           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 175           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 176           write (iout,*) "The internal geometry is:"
 177 c          call intout
 178 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 179           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 180           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 181           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 182           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 183           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 184           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 185           write (iout,*) "The components of the energy are:"
 186           call enerprint(energia(0),fT)
 187           write (iout,*)
 188      &      "This conformation WILL NOT be added to the database."
 189           call flush(iout)
 190           goto 121
 191         else
 192 #ifdef DEBUG
 193           if (ipar.eq.iparm) write (iout,*) i,iparm,
 194      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 195 #endif
 196           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 197      &      dabs(eini-energia(0)).gt.tole) then
 198             if (errmsg_count.le.maxerrmsg_count) then
 199               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 200      &         "Warning: energy differs remarkably from ",
 201      &         " the value read in: ",energia(0),eini," point",
 202      &         iii+1,indstart(me1)+iii," T",
 203      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 204               errmsg_count=errmsg_count+1
 205               if (errmsg_count.gt.maxerrmsg_count)
 206      &          write (iout,*) "Too many warning messages"
 207               if (einicheck.gt.1) then
 208                 write (iout,*) "Calculation stopped."
 209                 call flush(iout)
 210 #ifdef MPI
 211                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 212 #endif
 213                 call flush(iout)
 214                 return1
 215               endif
 216             endif
 217           endif
 218           potE(iii+1,iparm)=energia(0)
 219           do k=1,21
 220             enetb(k,iii+1,iparm)=energia(k)
 221           enddo
 222 #ifdef DEBUG
 223           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
 224      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 225 c          call enerprint(energia(0),fT)
 226 #endif
 227 #ifdef DEBUG
 228           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 229           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 230           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 231           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 232           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 233           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 234           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 235           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 236           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 237           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 238           write (iout,'(f10.5,i10)') rmsdev,iscor
 239           call enerprint(energia(0),fT)
 240         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 241 #endif
 242         endif
 243
 244         enddo ! iparm
 245
 246         iii=iii+1
 247         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 248         write (ientout,rec=iii)
 249      &   ((csingle(l,k),l=1,3),k=1,nres),
 250      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 251      &   nss,(ihpb(k),jhpb(k),k=1,nss),
 252      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 253 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 254 #ifdef MPI
 255         if (separate_parset) then
 256           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 257         else
 258           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 259         endif
 260 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 261 c     &   " snk",snk_p(iR,ib,ipar)
 262 #else
 263         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 264 #endif
 265   121   continue
 266       enddo
 267 #ifdef MPI
 268       scount(me)=iii
 269       write (iout,*) "Me",me," scount",scount(me)
 270       call flush(iout)
 271 c  Master gathers updated numbers of conformations written by all procs.
 272       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount(0), 1,
 273      &  MPI_INTEGER, WHAM_COMM, IERROR)
 274       indstart(0)=1
 275       indend(0)=scount(0)
 276       do i=1, Nprocs-1
 277         indstart(i)=indend(i-1)+1
 278         indend(i)=indstart(i)+scount(i)-1
 279       enddo
 280       write (iout,*)
 281       write (iout,*) "Revised conformation counts"
 282       do i=0,nprocs1-1
 283         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 284      &    "Processor",i," indstart",indstart(i),
 285      &    " indend",indend(i)," count",scount(i)
 286       enddo
 287       call flush(iout)
 288       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 289      &  MaxR*MaxT_h*nParmSet,
 290      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 291 #endif
 292       stot(islice)=0
 293       do iparm=1,nParmSet
 294         do ib=1,nT_h(iparm)
 295           do i=1,nR(ib,iparm)
 296             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 297           enddo
 298         enddo
 299       enddo
 300       write (iout,*) "Revised SNK"
 301       do iparm=1,nParmSet
 302         do ib=1,nT_h(iparm)
 303           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 304      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 305      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 306           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 307         enddo
 308       enddo
 309       write (iout,'("Total",i10)') stot(islice)
 310       call flush(iout)
 311       return
 312   101 write (iout,*) "Error in scratchfile."
 313       call flush(iout)
 314       return1
 315       end
 316 c------------------------------------------------------------------------------
 317       subroutine write_dbase(islice,*)
 318       implicit none
 319       include "DIMENSIONS"
 320       include "DIMENSIONS.ZSCOPT"
 321       include "DIMENSIONS.FREE"
 322       include "DIMENSIONS.COMPAR"
 323 #ifdef MPI
 324       include "mpif.h"
 325       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 326       include "COMMON.MPI"
 327 #endif
 328       include "COMMON.CONTROL"
 329       include "COMMON.CHAIN"
 330       include "COMMON.IOUNITS"
 331       include "COMMON.PROTFILES"
 332       include "COMMON.NAMES"
 333       include "COMMON.VAR"
 334       include "COMMON.SBRIDGE"
 335       include "COMMON.GEO"
 336       include "COMMON.FFIELD"
 337       include "COMMON.ENEPS"
 338       include "COMMON.LOCAL"
 339       include "COMMON.WEIGHTS"
 340       include "COMMON.INTERACT"
 341       include "COMMON.FREE"
 342       include "COMMON.ENERGIES"
 343       include "COMMON.COMPAR"
 344       include "COMMON.PROT"
 345       character*64 nazwa
 346       character*80 bxname,cxname
 347       character*64 bprotfile_temp
 348       character*3 liczba,licz
 349       character*2 licz2
 350       integer i,itj,ii,iii,j,k,l
 351       integer ixdrf,iret
 352       integer iscor,islice
 353       double precision rmsdev,efree,eini
 354       real*4 csingle(3,maxres2)
 355       double precision energ
 356       integer ilen,iroof
 357       external ilen,iroof
 358       integer ir,ib,iparm
 359       write (licz2,'(bz,i2.2)') islice
 360       call opentmp(islice,ientout,bprotfile_temp)
 361       write (iout,*) "bprotfile_temp ",bprotfile_temp
 362       call flush(iout)
 363       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 364      &   .and. ensembles.eq.0) then
 365         close(ientout,status="delete")
 366         return
 367       endif
 368 #ifdef MPI
 369       write (liczba,'(bz,i3.3)') me
 370       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 371         if (.not.separate_parset) then
 372           bxname = prefix(:ilen(prefix))//liczba//".bx"
 373         else
 374           write (licz,'(bz,i3.3)') myparm
 375           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 376         endif
 377         open (ientin,file=bxname,status="unknown",
 378      &    form="unformatted",access="direct",recl=lenrec1)
 379       endif
 380 #else
 381       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 382         if (nslice.eq.1) then
 383           bxname = prefix(:ilen(prefix))//".bx"
 384         else
 385           bxname = prefix(:ilen(prefix))//
 386      &           "_slice_"//licz2//".bx"
 387         endif
 388         open (ientin,file=bxname,status="unknown",
 389      &    form="unformatted",access="direct",recl=lenrec1)
 390         write (iout,*) "Calculating energies; writing geometry",
 391      & " and energy components to ",bxname(:ilen(bxname))
 392       endif
 393 #if (defined(AIX) && !defined(JUBL))
 394         call xdrfopen_(ixdrf,cxname, "w", iret)
 395 #else
 396         call xdrfopen(ixdrf,cxname, "w", iret)
 397 #endif
 398         if (iret.eq.0) then
 399           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 400           cxfile=.fale.
 401         endif
 402       endif
 403 #endif
 404       if (indpdb.gt.0) then
 405         if (nslice.eq.1) then
 406 #ifdef MPI
 407          if (.not.separate_parset) then
 408            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 409      &       //liczba//'.stat'
 410          else
 411            write (licz,'(bz,i3.3)') myparm
 412            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 413      &      pot(:ilen(pot))//liczba//'.stat'
 414          endif
 415
 416 #else
 417           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 418 #endif
 419         else
 420 #ifdef MPI
 421          if (.not.separate_parset) then
 422           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 423      &      "_slice_"//licz2//liczba//'.stat'
 424          else
 425           write (licz,'(bz,i3.3)') myparm
 426           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 427      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 428          endif
 429 #else
 430           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 431      &      //"_slice_"//licz2//'.stat'
 432 #endif
 433         endif
 434         open(istat,file=statname,status="unknown")
 435       endif
 436
 437 #ifdef MPI
 438       do i=1,scount(me)
 439 #else
 440       do i=1,ntot(islice)
 441 #endif
 442         read(ientout,rec=i,err=101)
 443      &    ((csingle(l,k),l=1,3),k=1,nres),
 444      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 445      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 446      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 447 c        write (iout,*) iR,ib,iparm,eini,efree
 448         do j=1,2*nres
 449           do k=1,3
 450             c(k,j)=csingle(k,j)
 451           enddo
 452         enddo
 453         call int_from_cart1(.false.)
 454         iscore=0
 455 c        write (iout,*) "Calling conf_compar",i
 456 c        call flush(iout)
 457         if (indpdb.gt.0) then
 458           call conf_compar(i,.false.,.true.)
 459         endif
 460 c        write (iout,*) "Exit conf_compar",i
 461 c        call flush(iout)
 462         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 463      &    ((csingle(l,k),l=1,3),k=1,nres),
 464      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 465      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 466 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 467      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 468 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 469 #ifndef MPI
 470         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 471      &    -entfac(i),rms_nat,iscore)
 472 #endif
 473       enddo
 474       close(ientout,status="delete")
 475       close(istat)
 476       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 477 #ifdef MPI
 478       call MPI_Barrier(WHAM_COMM,IERROR)
 479       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 480      &   .and. ensembles.eq.0) return
 481       write (iout,*)
 482       if (bxfile .or. ensembles.gt.0) then
 483         if (nslice.eq.1) then
 484           if (.not.separate_parset) then
 485             bxname = prefix(:ilen(prefix))//".bx"
 486           else
 487             write (licz,'(bz,i3.3)') myparm
 488             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 489           endif
 490         else
 491           if (.not.separate_parset) then
 492             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 493           else
 494             write (licz,'(bz,i3.3)') myparm
 495             bxname = prefix(:ilen(prefix))//"par_"//licz//
 496      &        "_slice_"//licz2//".bx"
 497           endif
 498         endif
 499         open (ientout,file=bxname,status="unknown",
 500      &      form="unformatted",access="direct",recl=lenrec1)
 501         write (iout,*) "Master is creating binary database ",
 502      &   bxname(:ilen(bxname))
 503       endif
 504       if (cxfile) then
 505         if (nslice.eq.1) then
 506           if (.not.separate_parset) then
 507             cxname = prefix(:ilen(prefix))//".cx"
 508           else
 509             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 510           endif
 511         else
 512           if (.not.separate_parset) then
 513             cxname = prefix(:ilen(prefix))//
 514      &             "_slice_"//licz2//".cx"
 515           else
 516             cxname = prefix(:ilen(prefix))//"_par"//licz//
 517      &             "_slice_"//licz2//".cx"
 518           endif
 519         endif
 520 #if (defined(AIX) && !defined(JUBL))
 521         call xdrfopen_(ixdrf,cxname, "w", iret)
 522 #else
 523         call xdrfopen(ixdrf,cxname, "w", iret)
 524 #endif
 525         if (iret.eq.0) then
 526           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 527           cxfile=.false.
 528         endif
 529       endif
 530       do j=0,nprocs-1
 531         write (liczba,'(bz,i3.3)') j
 532         if (separate_parset) then
 533           write (licz,'(bz,i3.3)') myparm
 534           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 535         else
 536           bxname = prefix(:ilen(prefix))//liczba//".bx"
 537         endif
 538         open (ientin,file=bxname,status="unknown",
 539      &    form="unformatted",access="direct",recl=lenrec1)
 540         write (iout,*) "Master is reading conformations from ",
 541      &   bxname(:ilen(bxname))
 542         iii = 0
 543 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 544 c        call flush(iout)
 545         do i=indstart(j),indend(j)
 546           iii = iii+1
 547           read(ientin,rec=iii,err=101)
 548      &      ((csingle(l,k),l=1,3),k=1,nres),
 549      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 550      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 551      &      eini,efree,rmsdev,iscor
 552           if (bxfile .or. ensembles.gt.0) then
 553             write (ientout,rec=i)
 554      &        ((csingle(l,k),l=1,3),k=1,nres),
 555      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 556      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 557      &        eini,efree,rmsdev,iscor
 558           endif
 559           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 560 #ifdef DEBUG
 561           do k=1,2*nres
 562             do l=1,3
 563               c(l,k)=csingle(l,k)
 564             enddo
 565           enddo
 566           call int_from_cart1(.false.)
 567           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 568           write (iout,*) "The Cartesian geometry is:"
 569           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 570           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 571           write (iout,*) "The internal geometry is:"
 572           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 573           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 574           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 575           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 576           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 577           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 578           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 579           write (iout,'(f10.5,i5)') rmsdev,iscor
 580 #endif
 581         enddo ! i
 582         write (iout,*) iii," conformations (from",indstart(j)," to",
 583      &   indend(j),") read from ",
 584      &   bxname(:ilen(bxname))
 585         close (ientin,status="delete")
 586       enddo ! j
 587       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 588 #if (defined(AIX) && !defined(JUBL))
 589       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 590 #else
 591       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 592 #endif
 593 #endif
 594       return
 595   101 write (iout,*) "Error in scratchfile."
 596       call flush(iout)
 597       return1
 598       end
 599 c-------------------------------------------------------------------------------
 600       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 601       implicit none
 602       include "DIMENSIONS"
 603       include "DIMENSIONS.ZSCOPT"
 604       include "DIMENSIONS.FREE"
 605       include "DIMENSIONS.COMPAR"
 606 #ifdef MPI
 607       include "mpif.h"
 608       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 609       include "COMMON.MPI"
 610 #endif
 611       include "COMMON.CONTROL"
 612       include "COMMON.CHAIN"
 613       include "COMMON.IOUNITS"
 614       include "COMMON.PROTFILES"
 615       include "COMMON.NAMES"
 616       include "COMMON.VAR"
 617       include "COMMON.SBRIDGE"
 618       include "COMMON.GEO"
 619       include "COMMON.FFIELD"
 620       include "COMMON.ENEPS"
 621       include "COMMON.LOCAL"
 622       include "COMMON.WEIGHTS"
 623       include "COMMON.INTERACT"
 624       include "COMMON.FREE"
 625       include "COMMON.ENERGIES"
 626       include "COMMON.COMPAR"
 627       include "COMMON.PROT"
 628       integer i,j,itmp,iscor,iret,ixdrf
 629       double precision rmsdev,efree,eini
 630       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 631       real*4 prec
 632
 633 c      write (iout,*) "cxwrite"
 634 c      call flush(iout)
 635       prec=10000.0
 636       do i=1,nres
 637        do j=1,3
 638         xoord(j,i)=csingle(j,i)
 639        enddo
 640       enddo
 641       do i=nnt,nct
 642        do j=1,3
 643         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 644        enddo
 645       enddo
 646
 647       itmp=nres+nct-nnt+1
 648
 649 c      write (iout,*) "itmp",itmp
 650 c      call flush(iout)
 651 #if (defined(AIX) && !defined(JUBL))
 652       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 653
 654 c      write (iout,*) "xdrf3dfcoord"
 655 c      call flush(iout)
 656       call xdrfint_(ixdrf, nss, iret)
 657       do j=1,nss
 658         call xdrfint_(ixdrf, ihpb(j), iret)
 659         call xdrfint_(ixdrf, jhpb(j), iret)
 660       enddo
 661       call xdrffloat_(ixdrf,real(eini),iret)
 662       call xdrffloat_(ixdrf,real(efree),iret)
 663       call xdrffloat_(ixdrf,real(rmsdev),iret)
 664       call xdrfint_(ixdrf,iscor,iret)
 665 #else
 666       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 667
 668       call xdrfint(ixdrf, nss, iret)
 669       do j=1,nss
 670         call xdrfint(ixdrf, ihpb(j), iret)
 671         call xdrfint(ixdrf, jhpb(j), iret)
 672       enddo
 673       call xdrffloat(ixdrf,real(eini),iret)
 674       call xdrffloat(ixdrf,real(efree),iret)
 675       call xdrffloat(ixdrf,real(rmsdev),iret)
 676       call xdrfint(ixdrf,iscor,iret)
 677 #endif
 678
 679       return
 680       end
 681 c------------------------------------------------------------------------------
 682       logical function conf_check(ii,iprint)
 683       implicit none
 684       include "DIMENSIONS"
 685       include "DIMENSIONS.ZSCOPT"
 686       include "DIMENSIONS.FREE"
 687 #ifdef MPI
 688       include "mpif.h"
 689       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 690       include "COMMON.MPI"
 691 #endif
 692       include "COMMON.CHAIN"
 693       include "COMMON.IOUNITS"
 694       include "COMMON.PROTFILES"
 695       include "COMMON.NAMES"
 696       include "COMMON.VAR"
 697       include "COMMON.SBRIDGE"
 698       include "COMMON.GEO"
 699       include "COMMON.FFIELD"
 700       include "COMMON.ENEPS"
 701       include "COMMON.LOCAL"
 702       include "COMMON.WEIGHTS"
 703       include "COMMON.INTERACT"
 704       include "COMMON.FREE"
 705       include "COMMON.ENERGIES"
 706       include "COMMON.CONTROL"
 707       include "COMMON.TORCNSTR"
 708       integer j,k,l,ii,itj,iprint
 709       if (.not.check_conf) then
 710         conf_check=.true.
 711         return
 712       endif
 713       call int_from_cart1(.false.)
 714       do j=nnt+1,nct
 715         if (itype(j-1).ne.21 .and. itype(j).ne.21 .and.
 716      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
 717           if (iprint.gt.0)
 718      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 719      &      " for conformation",ii
 720           if (iprint.gt.1) then
 721             write (iout,*) "The Cartesian geometry is:"
 722             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 723             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 724             write (iout,*) "The internal geometry is:"
 725             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 726             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 727             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 728             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 729             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 730             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 731           endif
 732           if (iprint.gt.0) write (iout,*)
 733      &      "This conformation WILL NOT be added to the database."
 734           conf_check=.false.
 735           return
 736         endif
 737       enddo
 738       do j=nnt,nct
 739         itj=itype(j)
 740         if (itype(j).ne.10 .and.itype(j).ne.21 .and.
 741      &     (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 742           if (iprint.gt.0)
 743      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 744      &     " for conformation",ii
 745           if (iprint.gt.1) then
 746             write (iout,*) "The Cartesian geometry is:"
 747             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 748             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 749             write (iout,*) "The internal geometry is:"
 750             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 751             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 752             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 753             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 754             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 755             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 756           endif
 757           if (iprint.gt.0) write (iout,*)
 758      &      "This conformation WILL NOT be added to the database."
 759           conf_check=.false.
 760           return
 761         endif
 762       enddo
 763       do j=3,nres
 764         if (theta(j).le.0.0d0) then
 765           if (iprint.gt.0)
 766      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 767           if (iprint.gt.1) then
 768             write (iout,*) "The Cartesian geometry is:"
 769             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 770             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 771             write (iout,*) "The internal geometry is:"
 772             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 773             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 774             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 775             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 776             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 777             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 778           endif
 779           if (iprint.gt.0) write (iout,*)
 780      &      "This conformation WILL NOT be added to the database."
 781           conf_check=.false.
 782           return
 783         endif
 784         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 785       enddo
 786       conf_check=.true.
 787 c      write (iout,*) "conf_check passed",ii
 788       return
 789       end