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