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