source/wham/src/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,scme,itmp
  37       integer ires
  38       integer ir,ib,ipar,iparm
  39       integer iscor,islice
  40       real*4 csingle(3,maxres2)
  41       double precision energ
  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 cc       write(iout,*), 'NAWEJ',i,eini
  80          if (indpdb.gt.0) then
  81            do k=1,nres
  82              do l=1,3
  83                c(l,k)=csingle(l,k)
  84              enddo
  85            enddo
  86            do k=nnt,nct
  87              do l=1,3
  88                c(l,k+nres)=csingle(l,k+nres)
  89              enddo
  90            enddo
  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         call enerprint(energia(0),fT)
 161 c        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 162 c        write (iout,*) "ftors",ftors
 163       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 164      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 165       do ires=1,nres
 166         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 167      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 168      &    (c(j,ires+nres),j=1,3)
 169       enddo
 170         call intout
 171 #endif
 172         if (energia(0).ge.1.0d20) then
 173           write (iout,*) "NaNs detected in some of the energy",
 174      &     " components for conformation",ii+1
 175           write (iout,*) "The Cartesian geometry is:"
 176           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 177           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 178           write (iout,*) "The internal geometry is:"
 179 c          call intout
 180 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 181           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 182           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 183           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 184           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 185           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 186           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 187           write (iout,*) "The components of the energy are:"
 188           call enerprint(energia(0),fT)
 189           write (iout,*)
 190      &      "This conformation WILL NOT be added to the database."
 191           call flush(iout)
 192           goto 121
 193         else
 194 #ifdef DEBUG
 195           if (ipar.eq.iparm) write (iout,*) i,iparm,
 196      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 197 #endif
 198           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 199      &      dabs(eini-energia(0)).gt.tole) then
 200             if (errmsg_count.le.maxerrmsg_count) then
 201               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 202      &         "Warning: energy differs remarkably from ",
 203      &         " the value read in: ",energia(0),eini," point",
 204      &         iii+1,indstart(me1)+iii," T",
 205      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 206              call enerprint(energia(0),fT)
 207             write (iout,'(4f10.5,2i5)') 0.0,energia(0),0.0,
 208      &       1.0d0/(beta_h(ib,ipar)*1.987D-3),
 209      &       0,0
 210             write(iout,'(8f10.5)')
 211      &       ((c(l,k),l=1,3),k=1,nres),
 212      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 213 c             itmp=ipdb
 214 c             ipdb=iout
 215 c             call pdbout(iii+1,beta_h(ib,ipar),
 216 c     &                   eini,energia(0),0.0d0,rmsdev)
 217              write (iout,*)
 218              ipdb=itmp
 219
 220               errmsg_count=errmsg_count+1
 221               if (errmsg_count.gt.maxerrmsg_count)
 222      &          write (iout,*) "Too many warning messages"
 223               if (einicheck.gt.1) then
 224                 write (iout,*) "Calculation stopped."
 225                 call flush(iout)
 226 #ifdef MPI
 227                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 228 #endif
 229                 call flush(iout)
 230                 return1
 231               endif
 232             endif
 233           endif
 234           potE(iii+1,iparm)=energia(0)
 235           do k=1,max_ene
 236             enetb(k,iii+1,iparm)=energia(k)
 237           enddo
 238 c          write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i)
 239 c          call enerprint(energia(0),fT)
 240 #ifdef DEBUG
 241           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
 242      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 243           call enerprint(energia(0),fT)
 244           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 245           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 246           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 247           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 248           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 249           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 250           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 251           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 252           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 253           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 254           write (iout,'(f10.5,i10)') rmsdev,iscor
 255           call enerprint(energia(0),fT)
 256         write(liczba,'(bz,i3.3)') me
 257         nazwa="test"//liczba//".pdb"
 258         write (iout,*) "pdb file",nazwa
 259         open (ipdb,file=nazwa,position="append")
 260         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 261         close(ipdb)
 262 #endif
 263         endif
 264
 265         enddo ! iparm
 266
 267         iii=iii+1
 268         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
 269         write (ientout,rec=iii)
 270      &   ((csingle(l,k),l=1,3),k=1,nres),
 271      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 272      &   nss,(ihpb(k),jhpb(k),k=1,nss),
 273      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 274 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 275 #ifdef MPI
 276         if (separate_parset) then
 277           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 278         else
 279           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 280         endif
 281 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 282 c     &   " snk",snk_p(iR,ib,ipar)
 283 #else
 284         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 285 #endif
 286   121   continue
 287       enddo
 288 #ifdef MPI
 289       scount(me)=iii
 290       write (iout,*) "Me",me," scount",scount(me)
 291       call flush(iout)
 292 c  Master gathers updated numbers of conformations written by all procs.
 293       scme = scount(me)
 294       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount(0), 1,
 295      &  MPI_INTEGER, WHAM_COMM, IERROR)
 296       indstart(0)=1
 297       indend(0)=scount(0)
 298       do i=1, Nprocs-1
 299         indstart(i)=indend(i-1)+1
 300         indend(i)=indstart(i)+scount(i)-1
 301       enddo
 302       write (iout,*)
 303       write (iout,*) "Revised conformation counts"
 304       do i=0,nprocs1-1
 305         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 306      &    "Processor",i," indstart",indstart(i),
 307      &    " indend",indend(i)," count",scount(i)
 308       enddo
 309       call flush(iout)
 310       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 311      &  MaxR*MaxT_h*nParmSet,
 312      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 313 #endif
 314       stot(islice)=0
 315       do iparm=1,nParmSet
 316         do ib=1,nT_h(iparm)
 317           do i=1,nR(ib,iparm)
 318             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 319           enddo
 320         enddo
 321       enddo
 322       write (iout,*) "Revised SNK"
 323       do iparm=1,nParmSet
 324         do ib=1,nT_h(iparm)
 325           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 326      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 327      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 328           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 329         enddo
 330       enddo
 331       write (iout,'("Total",i10)') stot(islice)
 332       call flush(iout)
 333       return
 334   101 write (iout,*) "Error in scratchfile."
 335       call flush(iout)
 336       return1
 337       end
 338 c------------------------------------------------------------------------------
 339       subroutine write_dbase(islice,*)
 340       implicit none
 341       include "DIMENSIONS"
 342       include "DIMENSIONS.ZSCOPT"
 343       include "DIMENSIONS.FREE"
 344       include "DIMENSIONS.COMPAR"
 345 #ifdef MPI
 346       include "mpif.h"
 347       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 348       include "COMMON.MPI"
 349 #endif
 350       include "COMMON.CONTROL"
 351       include "COMMON.CHAIN"
 352       include "COMMON.IOUNITS"
 353       include "COMMON.PROTFILES"
 354       include "COMMON.NAMES"
 355       include "COMMON.VAR"
 356       include "COMMON.SBRIDGE"
 357       include "COMMON.GEO"
 358       include "COMMON.FFIELD"
 359       include "COMMON.ENEPS"
 360       include "COMMON.LOCAL"
 361       include "COMMON.WEIGHTS"
 362       include "COMMON.INTERACT"
 363       include "COMMON.FREE"
 364       include "COMMON.ENERGIES"
 365       include "COMMON.COMPAR"
 366       include "COMMON.PROT"
 367       character*64 nazwa
 368       character*80 bxname,cxname
 369       character*64 bprotfile_temp
 370       character*3 liczba,licz
 371       character*2 licz2
 372       integer i,itj,ii,iii,j,k,l
 373       integer ixdrf,iret
 374       integer iscor,islice
 375       double precision rmsdev,efree,eini
 376       real*4 csingle(3,maxres2)
 377       double precision energ
 378       integer ilen,iroof
 379       external ilen,iroof
 380       integer ir,ib,iparm
 381       write (licz2,'(bz,i2.2)') islice
 382       call opentmp(islice,ientout,bprotfile_temp)
 383       write (iout,*) "bprotfile_temp ",bprotfile_temp
 384       call flush(iout)
 385       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 386      &   .and. ensembles.eq.0) then
 387         close(ientout,status="delete")
 388         return
 389       endif
 390 #ifdef MPI
 391       write (liczba,'(bz,i3.3)') me
 392       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 393         if (.not.separate_parset) then
 394           bxname = prefix(:ilen(prefix))//liczba//".bx"
 395         else
 396           write (licz,'(bz,i3.3)') myparm
 397           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 398         endif
 399         open (ientin,file=bxname,status="unknown",
 400      &    form="unformatted",access="direct",recl=lenrec1)
 401       endif
 402 #else
 403       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 404         if (nslice.eq.1) then
 405           bxname = prefix(:ilen(prefix))//".bx"
 406         else
 407           bxname = prefix(:ilen(prefix))//
 408      &           "_slice_"//licz2//".bx"
 409         endif
 410         open (ientin,file=bxname,status="unknown",
 411      &    form="unformatted",access="direct",recl=lenrec1)
 412         write (iout,*) "Calculating energies; writing geometry",
 413      & " and energy components to ",bxname(:ilen(bxname))
 414       endif
 415 #if (defined(AIX) && !defined(JUBL))
 416         call xdrfopen_(ixdrf,cxname, "w", iret)
 417 #else
 418         call xdrfopen(ixdrf,cxname, "w", iret)
 419 #endif
 420         if (iret.eq.0) then
 421           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 422           cxfile=.fale.
 423         endif
 424       endif
 425 #endif
 426       if (indpdb.gt.0) then
 427         if (nslice.eq.1) then
 428 #ifdef MPI
 429          if (.not.separate_parset) then
 430            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 431      &       //liczba//'.stat'
 432          else
 433            write (licz,'(bz,i3.3)') myparm
 434            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 435      &      pot(:ilen(pot))//liczba//'.stat'
 436          endif
 437
 438 #else
 439           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 440 #endif
 441         else
 442 #ifdef MPI
 443          if (.not.separate_parset) then
 444           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 445      &      "_slice_"//licz2//liczba//'.stat'
 446          else
 447           write (licz,'(bz,i3.3)') myparm
 448           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 449      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 450          endif
 451 #else
 452           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 453      &      //"_slice_"//licz2//'.stat'
 454 #endif
 455         endif
 456         open(istat,file=statname,status="unknown")
 457       endif
 458
 459 #ifdef MPI
 460       do i=1,scount(me)
 461 #else
 462       do i=1,ntot(islice)
 463 #endif
 464 cc        if (dyn_ss) then
 465 cc        read(ientout,rec=i,err=101)
 466 cc     &    ((csingle(l,k),l=1,3),k=1,nres),
 467 cc     &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 468 cc     &    nss,(idssb(k),jdssb(k),k=1,nss),
 469 cc     &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 470 cc        idssb(k)=idssb(k)-nres
 471 cc        jdssb(k)=jdssb(k)-nres
 472 cc        else
 473         read(ientout,rec=i,err=101)
 474      &    ((csingle(l,k),l=1,3),k=1,nres),
 475      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 476      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 477      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 478 cc         endif
 479 cc        write (iout,*) 'CC', iR,ib,iparm,eini,efree
 480         do j=1,2*nres
 481           do k=1,3
 482             c(k,j)=csingle(k,j)
 483           enddo
 484         enddo
 485         call int_from_cart1(.false.)
 486         iscore=0
 487         if (indpdb.gt.0) then
 488           call conf_compar(i,.false.,.true.)
 489         endif
 490 c        if (dyn_ss) then
 491         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 492      &    ((csingle(l,k),l=1,3),k=1,nres),
 493      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 494      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 495 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 496      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 497 c        else
 498           if (bxfile .or.cxfile .or. ensembles.gt.0) write
 499      &     (ientin,rec=i)
 500      &    ((csingle(l,k),l=1,3),k=1,nres),
 501      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 502      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 503 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 504      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 505 c        endif
 506         write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 507 #ifndef MPI
 508         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 509      &    -entfac(i),rms_nat,iscore)
 510 #endif
 511       enddo
 512       close(ientout,status="delete")
 513       close(istat)
 514       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 515 #ifdef MPI
 516       call MPI_Barrier(WHAM_COMM,IERROR)
 517       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 518      &   .and. ensembles.eq.0) return
 519       write (iout,*)
 520       if (bxfile .or. ensembles.gt.0) then
 521         if (nslice.eq.1) then
 522           if (.not.separate_parset) then
 523             bxname = prefix(:ilen(prefix))//".bx"
 524           else
 525             write (licz,'(bz,i3.3)') myparm
 526             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 527           endif
 528         else
 529           if (.not.separate_parset) then
 530             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 531           else
 532             write (licz,'(bz,i3.3)') myparm
 533             bxname = prefix(:ilen(prefix))//"par_"//licz//
 534      &        "_slice_"//licz2//".bx"
 535           endif
 536         endif
 537         open (ientout,file=bxname,status="unknown",
 538      &      form="unformatted",access="direct",recl=lenrec1)
 539         write (iout,*) "Master is creating binary database ",
 540      &   bxname(:ilen(bxname))
 541       endif
 542       if (cxfile) then
 543         if (nslice.eq.1) then
 544           if (.not.separate_parset) then
 545             cxname = prefix(:ilen(prefix))//".cx"
 546           else
 547             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 548           endif
 549         else
 550           if (.not.separate_parset) then
 551             cxname = prefix(:ilen(prefix))//
 552      &             "_slice_"//licz2//".cx"
 553           else
 554             cxname = prefix(:ilen(prefix))//"_par"//licz//
 555      &             "_slice_"//licz2//".cx"
 556           endif
 557         endif
 558 #if (defined(AIX) && !defined(JUBL))
 559         call xdrfopen_(ixdrf,cxname, "w", iret)
 560 #else
 561         call xdrfopen(ixdrf,cxname, "w", iret)
 562 #endif
 563         if (iret.eq.0) then
 564           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 565           cxfile=.false.
 566         endif
 567       endif
 568       do j=0,nprocs-1
 569         write (liczba,'(bz,i3.3)') j
 570         if (separate_parset) then
 571           write (licz,'(bz,i3.3)') myparm
 572           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 573         else
 574           bxname = prefix(:ilen(prefix))//liczba//".bx"
 575         endif
 576         open (ientin,file=bxname,status="unknown",
 577      &    form="unformatted",access="direct",recl=lenrec1)
 578         write (iout,*) "Master is reading conformations from ",
 579      &   bxname(:ilen(bxname))
 580         iii = 0
 581 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 582 c        call flush(iout)
 583         do i=indstart(j),indend(j)
 584           iii = iii+1
 585 cc          if (dyn_ss) then
 586 cc          read(ientin,rec=iii,err=101)
 587 cc     &      ((csingle(l,k),l=1,3),k=1,nres),
 588 cc     &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 589 cc     &      nss,(idssb(k),jdssb(k),k=1,nss),
 590 cc    &      eini,efree,rmsdev,iscor
 591 cc          idssb(k)=idssb(k)-nres
 592 cc          jdssb(k)=jdssb(k)-nres
 593 cc          else
 594           read(ientin,rec=iii,err=101)
 595      &      ((csingle(l,k),l=1,3),k=1,nres),
 596      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 597      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 598      &      eini,efree,rmsdev,iscor
 599 cc          endif
 600           if (bxfile .or. ensembles.gt.0) then
 601 cc          if (dyn_ss) then
 602 cc            write (ientout,rec=i)
 603 cc     &        ((csingle(l,k),l=1,3),k=1,nres),
 604 cc     &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 605 cc     &        nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss),
 606 cc    &        eini,efree,rmsdev,iscor
 607 cc           else
 608                         write (ientout,rec=i)
 609      &        ((csingle(l,k),l=1,3),k=1,nres),
 610      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 611      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 612      &        eini,efree,rmsdev,iscor
 613 cc           write(iout,*) "W poszukiwaniu zlotych galotow"
 614 cc           write(iout,*) "efree=",efree,iii
 615 cc           endif
 616           endif
 617           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 618 #ifdef DEBUG
 619           do k=1,2*nres
 620             do l=1,3
 621               c(l,k)=csingle(l,k)
 622             enddo
 623           enddo
 624           call int_from_cart1(.false.)
 625           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 626           write (iout,*) "The Cartesian geometry is:"
 627           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 628           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 629           write (iout,*) "The internal geometry is:"
 630           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 631           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 632           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 633           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 634           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 635           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 636           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 637           write (iout,'(f10.5,i5)') rmsdev,iscor
 638 #endif
 639         enddo ! i
 640         write (iout,*) iii," conformations (from",indstart(j)," to",
 641      &   indend(j),") read from ",
 642      &   bxname(:ilen(bxname))
 643         close (ientin,status="delete")
 644       enddo ! j
 645       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 646 #if (defined(AIX) && !defined(JUBL))
 647       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 648 #else
 649       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 650 #endif
 651 #endif
 652       return
 653   101 write (iout,*) "Error in scratchfile."
 654       call flush(iout)
 655       return1
 656       end
 657 c-------------------------------------------------------------------------------
 658       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 659       implicit none
 660       include "DIMENSIONS"
 661       include "DIMENSIONS.ZSCOPT"
 662       include "DIMENSIONS.FREE"
 663       include "DIMENSIONS.COMPAR"
 664 #ifdef MPI
 665       include "mpif.h"
 666       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 667       include "COMMON.MPI"
 668 #endif
 669       include "COMMON.CONTROL"
 670       include "COMMON.CHAIN"
 671       include "COMMON.IOUNITS"
 672       include "COMMON.PROTFILES"
 673       include "COMMON.NAMES"
 674       include "COMMON.VAR"
 675       include "COMMON.SBRIDGE"
 676       include "COMMON.GEO"
 677       include "COMMON.FFIELD"
 678       include "COMMON.ENEPS"
 679       include "COMMON.LOCAL"
 680       include "COMMON.WEIGHTS"
 681       include "COMMON.INTERACT"
 682       include "COMMON.FREE"
 683       include "COMMON.ENERGIES"
 684       include "COMMON.COMPAR"
 685       include "COMMON.PROT"
 686       integer i,j,itmp,iscor,iret,ixdrf
 687       double precision rmsdev,efree,eini
 688       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 689       real*4 prec
 690
 691 c      write (iout,*) "cxwrite"
 692 c      call flush(iout)
 693       prec=10000.0
 694       do i=1,nres
 695        do j=1,3
 696         xoord(j,i)=csingle(j,i)
 697        enddo
 698       enddo
 699       do i=nnt,nct
 700        do j=1,3
 701         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 702        enddo
 703       enddo
 704
 705       itmp=nres+nct-nnt+1
 706
 707 c      write (iout,*) "itmp",itmp
 708 c      call flush(iout)
 709 c       write (iout,*) "CNZ",eini,dyn_ss
 710 #if (defined(AIX) && !defined(JUBL))
 711       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 712
 713 c      write (iout,*) "xdrf3dfcoord"
 714 c      call flush(iout)
 715       call xdrfint_(ixdrf, nss, iret)
 716       do j=1,nss
 717 cc       if (dyn_ss) then
 718 cc       call xdrfint_(ixdrf, idssb(j)+nres, iret)
 719 cc        call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 720 cc       else
 721         call xdrfint_(ixdrf, ihpb(j), iret)
 722         call xdrfint_(ixdrf, jhpb(j), iret)
 723 cc       endif
 724       enddo
 725       call xdrffloat_(ixdrf,real(eini),iret)
 726       call xdrffloat_(ixdrf,real(efree),iret)
 727       call xdrffloat_(ixdrf,real(rmsdev),iret)
 728       call xdrfint_(ixdrf,iscor,iret)
 729 #else
 730       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 731
 732       call xdrfint(ixdrf, nss, iret)
 733       do j=1,nss
 734 cc       if (dyn_ss) then
 735 cc        call xdrfint(ixdrf, idssb(j), iret)
 736 cc        call xdrfint(ixdrf, jdssb(j), iret)
 737 cc        idssb(j)=idssb(j)-nres
 738 cc       jdssb(j)=jdssb(j)-nres
 739 cc       else
 740         call xdrfint(ixdrf, ihpb(j), iret)
 741         call xdrfint(ixdrf, jhpb(j), iret)
 742 cc       endif
 743       enddo
 744       call xdrffloat(ixdrf,real(eini),iret)
 745       call xdrffloat(ixdrf,real(efree),iret)
 746       call xdrffloat(ixdrf,real(rmsdev),iret)
 747       call xdrfint(ixdrf,iscor,iret)
 748 #endif
 749
 750       return
 751       end
 752 c------------------------------------------------------------------------------
 753       logical function conf_check(ii,iprint)
 754       implicit none
 755       include "DIMENSIONS"
 756       include "DIMENSIONS.ZSCOPT"
 757       include "DIMENSIONS.FREE"
 758 #ifdef MPI
 759       include "mpif.h"
 760       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 761       include "COMMON.MPI"
 762 #endif
 763       include "COMMON.CHAIN"
 764       include "COMMON.IOUNITS"
 765       include "COMMON.PROTFILES"
 766       include "COMMON.NAMES"
 767       include "COMMON.VAR"
 768       include "COMMON.SBRIDGE"
 769       include "COMMON.GEO"
 770       include "COMMON.FFIELD"
 771       include "COMMON.ENEPS"
 772       include "COMMON.LOCAL"
 773       include "COMMON.WEIGHTS"
 774       include "COMMON.INTERACT"
 775       include "COMMON.FREE"
 776       include "COMMON.ENERGIES"
 777       include "COMMON.CONTROL"
 778       include "COMMON.TORCNSTR"
 779       integer j,k,l,ii,itj,iprint
 780       if (.not.check_conf) then
 781         conf_check=.true.
 782         return
 783       endif
 784       call int_from_cart1(.false.)
 785       do j=nnt+1,nct
 786         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
 787           if (iprint.gt.0)
 788      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 789      &      " for conformation",ii
 790           if (iprint.gt.1) then
 791             write (iout,*) "The Cartesian geometry is:"
 792             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 793             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 794             write (iout,*) "The internal geometry is:"
 795             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 796             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 797             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 798             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 799             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 800             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 801           endif
 802           if (iprint.gt.0) write (iout,*)
 803      &      "This conformation WILL NOT be added to the database."
 804           conf_check=.false.
 805           return
 806         endif
 807       enddo
 808       do j=nnt,nct
 809         itj=itype(j)
 810         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
 811           if (iprint.gt.0)
 812      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 813      &     " for conformation",ii
 814           if (iprint.gt.1) then
 815             write (iout,*) "The Cartesian geometry is:"
 816             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 817             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 818             write (iout,*) "The internal geometry is:"
 819             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 820             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 821             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 822             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 823             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 824             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 825           endif
 826           if (iprint.gt.0) write (iout,*)
 827      &      "This conformation WILL NOT be added to the database."
 828           conf_check=.false.
 829           return
 830         endif
 831       enddo
 832       do j=3,nres
 833         if (theta(j).le.0.0d0) then
 834           if (iprint.gt.0)
 835      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 836           if (iprint.gt.1) then
 837             write (iout,*) "The Cartesian geometry is:"
 838             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 839             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 840             write (iout,*) "The internal geometry is:"
 841             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 842             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 843             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 844             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 845             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 846             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 847           endif
 848           if (iprint.gt.0) write (iout,*)
 849      &      "This conformation WILL NOT be added to the database."
 850           conf_check=.false.
 851           return
 852         endif
 853         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 854       enddo
 855       conf_check=.true.
 856 c      write (iout,*) "conf_check passed",ii
 857       return
 858       end