source/unres/src_MD-M-SAXS.safe/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read job setup parameters
  14       call read_control
  15 C Read force-field parameters except weights
  16       call parmread
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 C Read molecule information, molecule geometry, energy-term weights, and
  38 C restraints if requested
  39       call molread
  40 C Print restraint information
  41 #ifdef MPI
  42       if (.not. out1file .or. me.eq.king) then
  43 #endif
  44       if (nhpb.gt.nss)  then
  45         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  46      &  "The following",nhpb-nss,
  47      &  " distance restraints have been imposed:",
  48      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  49      &  "  score"," type"
  50         do i=nss+1,nhpb
  51           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  52      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  53      &  irestr_type(i)
  54         enddo
  55       endif
  56       if (npeak.gt.0) then
  57         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  58      &  "The following",npeak,
  59      &  " NMR peak restraints have been imposed:",
  60      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  61      &  "  score"," type"," ipeak"
  62         do i=1,npeak
  63           do j=ipeak(1,i),ipeak(2,i)
  64             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  65      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  66      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  67           enddo
  68         enddo
  69         write (iout,*) "The ipeak array"
  70         do i=1,npeak
  71           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  72         enddo
  73       endif
  74 #ifdef MPI
  75       endif
  76 #endif
  77 c      print *,"Processor",myrank," leaves READRTNS"
  78       return
  79       end
  80 C-------------------------------------------------------------------------------
  81       subroutine read_control
  82 C
  83 C Read contorl data
  84 C
  85       implicit real*8 (a-h,o-z)
  86       include 'DIMENSIONS'
  87 #ifdef MP
  88       include 'mpif.h'
  89       logical OKRandom, prng_restart
  90       real*8  r1
  91 #endif
  92       include 'COMMON.IOUNITS'
  93       include 'COMMON.TIME1'
  94       include 'COMMON.THREAD'
  95       include 'COMMON.SBRIDGE'
  96       include 'COMMON.CONTROL'
  97       include 'COMMON.MCM'
  98       include 'COMMON.MAP'
  99       include 'COMMON.HEADER'
 100       include 'COMMON.CSA'
 101       include 'COMMON.CHAIN'
 102       include 'COMMON.MUCA'
 103       include 'COMMON.MD'
 104       include 'COMMON.FFIELD'
 105       include 'COMMON.INTERACT'
 106       include 'COMMON.SETUP'
 107       include 'COMMON.SPLITELE'
 108       include 'COMMON.SHIELD'
 109       include 'COMMON.GEO'
 110       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 111       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 112       character*80 ucase
 113       character*320 controlcard
 114
 115       nglob_csa=0
 116       eglob_csa=1d99
 117       nmin_csa=0
 118       read (INP,'(a)') titel
 119       call card_concat(controlcard)
 120 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 121 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 122       call reada(controlcard,'SEED',seed,0.0D0)
 123       call random_init(seed)
 124 C Set up the time limit (caution! The time must be input in minutes!)
 125       read_cart=index(controlcard,'READ_CART').gt.0
 126       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 127 C this variable with_theta_constr is the variable which allow to read and execute the
 128 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 129       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 130       write (iout,*) "with_theta_constr ",with_theta_constr
 131       call readi(controlcard,'NSAXS',nsaxs,0)
 132       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 133       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 134       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 135       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 136      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 137       call readi(controlcard,'SYM',symetr,1)
 138       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 139       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 140       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 141       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 142       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 143       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 144       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 145       call reada(controlcard,'DRMS',drms,0.1D0)
 146       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 147        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 148        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 149        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 150        write (iout,'(a,f10.1)')'DRMS    = ',drms
 151        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 152        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 153       endif
 154       call readi(controlcard,'NZ_START',nz_start,0)
 155       call readi(controlcard,'NZ_END',nz_end,0)
 156       call readi(controlcard,'IZ_SC',iz_sc,0)
 157       timlim=60.0D0*timlim
 158       safety = 60.0d0*safety
 159       timem=timlim
 160       modecalc=0
 161       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 162       minim=(index(controlcard,'MINIMIZE').gt.0)
 163       dccart=(index(controlcard,'CART').gt.0)
 164       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 165       overlapsc=.not.overlapsc
 166       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 167       searchsc=.not.searchsc
 168       sideadd=(index(controlcard,'SIDEADD').gt.0)
 169       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 170       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 171       outpdb=(index(controlcard,'PDBOUT').gt.0)
 172       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 173       pdbref=(index(controlcard,'PDBREF').gt.0)
 174       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 175       indpdb=index(controlcard,'PDBSTART')
 176       extconf=(index(controlcard,'EXTCONF').gt.0)
 177       AFMlog=(index(controlcard,'AFM'))
 178       selfguide=(index(controlcard,'SELFGUIDE'))
 179 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 180       call readi(controlcard,'TUBEMOD',tubelog,0)
 181 c      write (iout,*) TUBElog,"TUBEMODE"
 182       call readi(controlcard,'IPRINT',iprint,0)
 183 C SHIELD keyword sets if the shielding effect of side-chains is used
 184 C 0 denots no shielding is used all peptide are equally despite the
 185 C solvent accesible area
 186 C 1 the newly introduced function
 187 C 2 reseved for further possible developement
 188       call readi(controlcard,'SHIELD',shield_mode,0)
 189 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 190         write(iout,*) "shield_mode",shield_mode
 191 C      endif
 192       call readi(controlcard,'TORMODE',tor_mode,0)
 193 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 194         write(iout,*) "torsional and valence angle mode",tor_mode
 195       call readi(controlcard,'MAXGEN',maxgen,10000)
 196       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 197       call readi(controlcard,"KDIAG",kdiag,0)
 198       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 199       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 200      & write (iout,*) "RESCALE_MODE",rescale_mode
 201       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 202       if (index(controlcard,'REGULAR').gt.0.0D0) then
 203         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 204         modecalc=1
 205         refstr=.true.
 206       endif
 207       if (index(controlcard,'CHECKGRAD').gt.0) then
 208         modecalc=5
 209         if (index(controlcard,' CART').gt.0) then
 210           icheckgrad=1
 211         elseif (index(controlcard,'CARINT').gt.0) then
 212           icheckgrad=2
 213         else
 214           icheckgrad=3
 215         endif
 216         call reada(controlcard,'DELTA',aincr,1.0d-7)
 217 c        write (iout,*) "icheckgrad",icheckgrad
 218       elseif (index(controlcard,'THREAD').gt.0) then
 219         modecalc=2
 220         call readi(controlcard,'THREAD',nthread,0)
 221         if (nthread.gt.0) then
 222           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 223         else
 224           if (fg_rank.eq.0)
 225      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 226           stop 'Error termination in Read_Control.'
 227         endif
 228       else if (index(controlcard,'MCMA').gt.0) then
 229         modecalc=3
 230       else if (index(controlcard,'MCEE').gt.0) then
 231         modecalc=6
 232       else if (index(controlcard,'MULTCONF').gt.0) then
 233         modecalc=4
 234       else if (index(controlcard,'MAP').gt.0) then
 235         modecalc=7
 236         call readi(controlcard,'MAP',nmap,0)
 237       else if (index(controlcard,'CSA').gt.0) then
 238         modecalc=8
 239 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 240 crc
 241 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 242 crc
 243 crc        modecalc=9
 244 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 245 cfmc        modecalc=10
 246       else if (index(controlcard,'SOFTREG').gt.0) then
 247         modecalc=11
 248       else if (index(controlcard,'CHECK_BOND').gt.0) then
 249         modecalc=-1
 250       else if (index(controlcard,'TEST').gt.0) then
 251         modecalc=-2
 252       else if (index(controlcard,'MD').gt.0) then
 253         modecalc=12
 254       else if (index(controlcard,'RE ').gt.0) then
 255         modecalc=14
 256       endif
 257
 258       lmuca=index(controlcard,'MUCA').gt.0
 259       call readi(controlcard,'MUCADYN',mucadyn,0)
 260       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 261       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 262      & then
 263        write (iout,*) 'MUCADYN=',mucadyn
 264        write (iout,*) 'MUCASMOOTH=',muca_smooth
 265       endif
 266
 267       iscode=index(controlcard,'ONE_LETTER')
 268       indphi=index(controlcard,'PHI')
 269       indback=index(controlcard,'BACK')
 270       iranconf=index(controlcard,'RAND_CONF')
 271       i2ndstr=index(controlcard,'USE_SEC_PRED')
 272       gradout=index(controlcard,'GRADOUT').gt.0
 273       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 274 C  DISTCHAINMAX become obsolete for periodic boundry condition
 275       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 276 C Reading the dimensions of box in x,y,z coordinates
 277       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 278       call reada(controlcard,'BOXY',boxysize,100.0d0)
 279       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 280 c Cutoff range for interactions
 281       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 282       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 283       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 284       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 285       if (lipthick.gt.0.0d0) then
 286        bordliptop=(boxzsize+lipthick)/2.0
 287        bordlipbot=bordliptop-lipthick
 288 C      endif
 289       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 290      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 291       buflipbot=bordlipbot+lipbufthick
 292       bufliptop=bordliptop-lipbufthick
 293       if ((lipbufthick*2.0d0).gt.lipthick)
 294      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 295       endif
 296       write(iout,*) "bordliptop=",bordliptop
 297       write(iout,*) "bordlipbot=",bordlipbot
 298       write(iout,*) "bufliptop=",bufliptop
 299       write(iout,*) "buflipbot=",buflipbot
 300       write (iout,*) "SHIELD MODE",shield_mode
 301       if (TUBElog.gt.0) then
 302        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 303        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 304        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 305        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 306        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 307        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 308        buftubebot=bordtubebot+tubebufthick
 309        buftubetop=bordtubetop-tubebufthick
 310       endif
 311 c      if (shield_mode.gt.0) then
 312 c      pi=3.141592d0
 313 C VSolvSphere the volume of solving sphere
 314 C      print *,pi,"pi"
 315 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 316 C there will be no distinction between proline peptide group and normal peptide
 317 C group in case of shielding parameters
 318 c      write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
 319 c      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 320 c      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 321 c      write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
 322 c     &  VSolvSphere_div
 323 C long axis of side chain
 324 c      do i=1,ntyp
 325 c      long_r_sidechain(i)=vbldsc0(1,i)
 326 c      short_r_sidechain(i)=sigma0(i)
 327 c      enddo
 328 c      buff_shield=1.0d0
 329 c      endif
 330       if (me.eq.king .or. .not.out1file )
 331      &  write (iout,*) "DISTCHAINMAX",distchainmax
 332
 333       if(me.eq.king.or..not.out1file)
 334      & write (iout,'(2a)') diagmeth(kdiag),
 335      &  ' routine used to diagonalize matrices.'
 336       return
 337       end
 338 c--------------------------------------------------------------------------
 339       subroutine read_REMDpar
 340 C
 341 C Read REMD settings
 342 C
 343       implicit real*8 (a-h,o-z)
 344       include 'DIMENSIONS'
 345       include 'COMMON.IOUNITS'
 346       include 'COMMON.TIME1'
 347       include 'COMMON.MD'
 348 #ifndef LANG0
 349       include 'COMMON.LANGEVIN'
 350 #else
 351       include 'COMMON.LANGEVIN.lang0'
 352 #endif
 353       include 'COMMON.INTERACT'
 354       include 'COMMON.NAMES'
 355       include 'COMMON.GEO'
 356       include 'COMMON.REMD'
 357       include 'COMMON.CONTROL'
 358       include 'COMMON.SETUP'
 359       character*80 ucase
 360       character*320 controlcard
 361       character*3200 controlcard1
 362       integer iremd_m_total
 363
 364       if(me.eq.king.or..not.out1file)
 365      & write (iout,*) "REMD setup"
 366
 367       call card_concat(controlcard)
 368       call readi(controlcard,"NREP",nrep,3)
 369       call readi(controlcard,"NSTEX",nstex,1000)
 370       call reada(controlcard,"RETMIN",retmin,10.0d0)
 371       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 372       mremdsync=(index(controlcard,'SYNC').gt.0)
 373       call readi(controlcard,"NSYN",i_sync_step,100)
 374       restart1file=(index(controlcard,'REST1FILE').gt.0)
 375       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 376       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 377       if(max_cache_traj_use.gt.max_cache_traj)
 378      &           max_cache_traj_use=max_cache_traj
 379       if(me.eq.king.or..not.out1file) then
 380 cd       if (traj1file) then
 381 crc caching is in testing - NTWX is not ignored
 382 cd        write (iout,*) "NTWX value is ignored"
 383 cd        write (iout,*) "  trajectory is stored to one file by master"
 384 cd        write (iout,*) "  before exchange at NSTEX intervals"
 385 cd       endif
 386        write (iout,*) "NREP= ",nrep
 387        write (iout,*) "NSTEX= ",nstex
 388        write (iout,*) "SYNC= ",mremdsync
 389        write (iout,*) "NSYN= ",i_sync_step
 390        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 391       endif
 392       remd_tlist=.false.
 393       if (index(controlcard,'TLIST').gt.0) then
 394          remd_tlist=.true.
 395          call card_concat(controlcard1)
 396          read(controlcard1,*) (remd_t(i),i=1,nrep)
 397          if(me.eq.king.or..not.out1file)
 398      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 399       endif
 400       remd_mlist=.false.
 401       if (index(controlcard,'MLIST').gt.0) then
 402          remd_mlist=.true.
 403          call card_concat(controlcard1)
 404          read(controlcard1,*) (remd_m(i),i=1,nrep)
 405          if(me.eq.king.or..not.out1file) then
 406           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 407           iremd_m_total=0
 408           do i=1,nrep
 409            iremd_m_total=iremd_m_total+remd_m(i)
 410           enddo
 411            write (iout,*) 'Total number of replicas ',iremd_m_total
 412           endif
 413          endif
 414         if (adaptive) then
 415           write (iout,*)
 416      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 417           call PMFread
 418         endif
 419       if(me.eq.king.or..not.out1file)
 420      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 421       return
 422       end
 423 c--------------------------------------------------------------------------
 424       subroutine read_MDpar
 425 C
 426 C Read MD settings
 427 C
 428       implicit real*8 (a-h,o-z)
 429       include 'DIMENSIONS'
 430       include 'COMMON.IOUNITS'
 431       include 'COMMON.TIME1'
 432       include 'COMMON.MD'
 433 #ifndef LANG0
 434       include 'COMMON.LANGEVIN'
 435 #else
 436       include 'COMMON.LANGEVIN.lang0'
 437 #endif
 438       include 'COMMON.INTERACT'
 439       include 'COMMON.NAMES'
 440       include 'COMMON.GEO'
 441       include 'COMMON.SETUP'
 442       include 'COMMON.CONTROL'
 443       include 'COMMON.SPLITELE'
 444       include 'COMMON.FFIELD'
 445       character*80 ucase
 446       character*320 controlcard
 447
 448       call card_concat(controlcard)
 449       call readi(controlcard,"NSTEP",n_timestep,1000000)
 450       call readi(controlcard,"NTWE",ntwe,100)
 451       call readi(controlcard,"NTWX",ntwx,1000)
 452       call reada(controlcard,"DT",d_time,1.0d-1)
 453       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 454       call reada(controlcard,"DAMAX",damax,1.0d1)
 455       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 456       call readi(controlcard,"LANG",lang,0)
 457       RESPA = index(controlcard,"RESPA") .gt. 0
 458       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 459       ntime_split0=ntime_split
 460       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 461       ntime_split0=ntime_split
 462 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 463 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 464       rest = index(controlcard,"REST").gt.0
 465       tbf = index(controlcard,"TBF").gt.0
 466       usampl = index(controlcard,"USAMPL").gt.0
 467       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 468       adaptive = index(controlcard,"ADAPTIVE").gt.0
 469       mdpdb = index(controlcard,"MDPDB").gt.0
 470       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 471       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 472       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 473       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 474       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 475       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 476       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 477       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 478       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 479       large = index(controlcard,"LARGE").gt.0
 480       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 481       rattle = index(controlcard,"RATTLE").gt.0
 482 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 483       nset=0
 484       if(usampl) then
 485         write (iout,*) "Umbrella sampling will be run"
 486         if (scale_umb.and.adaptive) then
 487           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 488           write (iout,*) "Select one of those and re-run the job."
 489           stop
 490         endif
 491         if (scale_umb) write (iout,*)
 492      &"Umbrella-restraint force constants will be scaled by temperature"
 493         call read_fragments
 494       endif
 495
 496       if(me.eq.king.or..not.out1file) then
 497        write (iout,*)
 498        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 499        write (iout,*)
 500        write (iout,'(a)') "The units are:"
 501        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 502        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 503      &  " acceleration: angstrom/(48.9 fs)**2"
 504        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 505        write (iout,*)
 506        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 507        write (iout,'(a60,f10.5,a)')
 508      &  "Initial time step of numerical integration:",d_time,
 509      &  " natural units"
 510        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 511        if (RESPA) then
 512         write (iout,'(2a,i4,a)')
 513      &    "A-MTS algorithm used; initial time step for fast-varying",
 514      &    " short-range forces split into",ntime_split," steps."
 515         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 516      &   r_cut," lambda",rlamb
 517        endif
 518        write (iout,'(2a,f10.5)')
 519      &  "Maximum acceleration threshold to reduce the time step",
 520      &  "/increase split number:",damax
 521        write (iout,'(2a,f10.5)')
 522      &  "Maximum predicted energy drift to reduce the timestep",
 523      &  "/increase split number:",edriftmax
 524        write (iout,'(a60,f10.5)')
 525      & "Maximum velocity threshold to reduce velocities:",dvmax
 526        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 527        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 528        if (rattle) write (iout,'(a60)')
 529      &  "Rattle algorithm used to constrain the virtual bonds"
 530       endif
 531       reset_fricmat=1000
 532       if (lang.gt.0) then
 533         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 534         call reada(controlcard,"RWAT",rwat,1.4d0)
 535         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 536         surfarea=index(controlcard,"SURFAREA").gt.0
 537         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 538         if(me.eq.king.or..not.out1file)then
 539          write (iout,'(/a,$)') "Langevin dynamics calculation"
 540          if (lang.eq.1) then
 541           write (iout,'(a/)')
 542      &      " with direct integration of Langevin equations"
 543          else if (lang.eq.2) then
 544           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 545          else if (lang.eq.3) then
 546           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 547          else if (lang.eq.4) then
 548           write (iout,'(a/)') " in overdamped mode"
 549          else
 550           write (iout,'(//a,i5)')
 551      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 552           stop
 553          endif
 554          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 555          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 556          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 557          write (iout,'(a60,f10.5)')
 558      &   "Scaling factor of the friction forces:",scal_fric
 559          if (surfarea) write (iout,'(2a,i10,a)')
 560      &     "Friction coefficients will be scaled by solvent-accessible",
 561      &     " surface area every",reset_fricmat," steps."
 562         endif
 563 c Calculate friction coefficients and bounds of stochastic forces
 564         eta=6*pi*cPoise*etawat
 565         if(me.eq.king.or..not.out1file)
 566      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 567      &   ,eta
 568         gamp=scal_fric*(pstok+rwat)*eta
 569         stdfp=dsqrt(2*Rb*t_bath/d_time)
 570         do i=1,ntyp
 571           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 572           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 573         enddo
 574         if(me.eq.king.or..not.out1file)then
 575          write (iout,'(/2a/)')
 576      &   "Radii of site types and friction coefficients and std's of",
 577      &   " stochastic forces of fully exposed sites"
 578          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 579          do i=1,ntyp
 580           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 581      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 582          enddo
 583         endif
 584       else if (tbf) then
 585         if(me.eq.king.or..not.out1file)then
 586          write (iout,'(a)') "Berendsen bath calculation"
 587          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 588          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 589          if (reset_moment)
 590      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 591      &   count_reset_moment," steps"
 592          if (reset_vel)
 593      &    write (iout,'(a,i10,a)')
 594      &    "Velocities will be reset at random every",count_reset_vel,
 595      &   " steps"
 596         endif
 597       else
 598         if(me.eq.king.or..not.out1file)
 599      &   write (iout,'(a31)') "Microcanonical mode calculation"
 600       endif
 601       if(me.eq.king.or..not.out1file)then
 602        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 603        if (usampl) then
 604           write(iout,*) "MD running with constraints."
 605           write(iout,*) "Equilibration time ", eq_time, " mtus."
 606           write(iout,*) "Constraining ", nfrag," fragments."
 607           write(iout,*) "Length of each fragment, weight and q0:"
 608           do iset=1,nset
 609            write (iout,*) "Set of restraints #",iset
 610            do i=1,nfrag
 611               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 612      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 613            enddo
 614            write(iout,*) "constraints between ", npair, "fragments."
 615            write(iout,*) "constraint pairs, weights and q0:"
 616            do i=1,npair
 617             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 618      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 619            enddo
 620            write(iout,*) "angle constraints within ", nfrag_back,
 621      &      "backbone fragments."
 622            if (loc_qlike) then
 623            write(iout,*) "fragment, weights, q0:"
 624            do i=1,nfrag_back
 625             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 626      &         ifrag_back(2,i,iset),
 627      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 628      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 629      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 630            enddo
 631            else
 632            write(iout,*) "fragment, weights:"
 633            do i=1,nfrag_back
 634             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 635      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 636      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 637            enddo
 638            endif
 639           enddo
 640         iset=mod(kolor,nset)+1
 641        endif
 642       endif
 643       if(me.eq.king.or..not.out1file)
 644      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 645       return
 646       end
 647 c------------------------------------------------------------------------------
 648       subroutine molread
 649 C
 650 C Read molecular data.
 651 C
 652       implicit real*8 (a-h,o-z)
 653       include 'DIMENSIONS'
 654 #ifdef MPI
 655       include 'mpif.h'
 656       integer error_msg
 657 #endif
 658       include 'COMMON.IOUNITS'
 659       include 'COMMON.GEO'
 660       include 'COMMON.VAR'
 661       include 'COMMON.INTERACT'
 662       include 'COMMON.LOCAL'
 663       include 'COMMON.NAMES'
 664       include 'COMMON.CHAIN'
 665       include 'COMMON.FFIELD'
 666       include 'COMMON.SBRIDGE'
 667       include 'COMMON.HEADER'
 668       include 'COMMON.CONTROL'
 669       include 'COMMON.DBASE'
 670       include 'COMMON.THREAD'
 671       include 'COMMON.CONTACTS'
 672       include 'COMMON.TORCNSTR'
 673       include 'COMMON.TIME1'
 674       include 'COMMON.BOUNDS'
 675       include 'COMMON.MD'
 676       include 'COMMON.SETUP'
 677       include 'COMMON.SHIELD'
 678       character*4 sequence(maxres)
 679       integer rescode
 680       double precision x(maxvar)
 681       character*256 pdbfile
 682       character*400 weightcard
 683       character*80 weightcard_t,ucase
 684       dimension itype_pdb(maxres)
 685       common /pizda/ itype_pdb
 686       logical seq_comp,fail
 687       double precision energia(0:n_ene)
 688       double precision secprob(3,maxdih_constr)
 689       integer ilen
 690       external ilen
 691 C
 692 C Body
 693 C
 694 C Read weights of the subsequent energy terms.
 695        call card_concat(weightcard)
 696        call reada(weightcard,'WLONG',wlong,1.0D0)
 697        call reada(weightcard,'WSC',wsc,wlong)
 698        call reada(weightcard,'WSCP',wscp,wlong)
 699        call reada(weightcard,'WELEC',welec,1.0D0)
 700        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 701        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 702        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 703        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 704        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 705        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 706        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 707        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 708        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 709        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 710        call reada(weightcard,'WBOND',wbond,1.0D0)
 711        call reada(weightcard,'WTOR',wtor,1.0D0)
 712        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 713        call reada(weightcard,'WANG',wang,1.0D0)
 714        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 715        call reada(weightcard,'SCAL14',scal14,0.4D0)
 716        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 717        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 718        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 719        call reada(weightcard,'TEMP0',temp0,300.0d0)
 720        call reada(weightcard,'WSHIELD',wshield,1.0d0)
 721        call reada(weightcard,'WLT',wliptran,0.0D0)
 722        call reada(weightcard,'WTUBE',wtube,1.0D0)
 723        call reada(weightcard,'WSAXS',wsaxs,1.0D0)
 724        if (index(weightcard,'SOFT').gt.0) ipot=6
 725 C 12/1/95 Added weight for the multi-body term WCORR
 726        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 727        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 728        weights(1)=wsc
 729        weights(2)=wscp
 730        weights(3)=welec
 731        weights(4)=wcorr
 732        weights(5)=wcorr5
 733        weights(6)=wcorr6
 734        weights(7)=wel_loc
 735        weights(8)=wturn3
 736        weights(9)=wturn4
 737        weights(10)=wturn6
 738        weights(11)=wang
 739        weights(12)=wscloc
 740        weights(13)=wtor
 741        weights(14)=wtor_d
 742        weights(15)=wstrain
 743        weights(16)=wvdwpp
 744        weights(17)=wbond
 745        weights(18)=scal14
 746        weights(21)=wsccor
 747        weights(26)=wsaxs
 748       if(me.eq.king.or..not.out1file)
 749      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 750      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 751      &  wturn4,wturn6
 752    10 format (/'Energy-term weights (unscaled):'//
 753      & 'WSCC=   ',f10.6,' (SC-SC)'/
 754      & 'WSCP=   ',f10.6,' (SC-p)'/
 755      & 'WELEC=  ',f10.6,' (p-p electr)'/
 756      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 757      & 'WBOND=  ',f10.6,' (stretching)'/
 758      & 'WANG=   ',f10.6,' (bending)'/
 759      & 'WSCLOC= ',f10.6,' (SC local)'/
 760      & 'WTOR=   ',f10.6,' (torsional)'/
 761      & 'WTORD=  ',f10.6,' (double torsional)'/
 762      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 763      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 764      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 765      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 766      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 767      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 768      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 769      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 770      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 771       if(me.eq.king.or..not.out1file)then
 772        if (wcorr4.gt.0.0d0) then
 773         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 774      &   'between contact pairs of peptide groups'
 775         write (iout,'(2(a,f5.3/))')
 776      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 777      &  'Range of quenching the correlation terms:',2*delt_corr
 778        else if (wcorr.gt.0.0d0) then
 779         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 780      &   'between contact pairs of peptide groups'
 781        endif
 782        write (iout,'(a,f8.3)')
 783      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 784        write (iout,'(a,f8.3)')
 785      &  'General scaling factor of SC-p interactions:',scalscp
 786       endif
 787       r0_corr=cutoff_corr-delt_corr
 788       do i=1,ntyp
 789         aad(i,1)=scalscp*aad(i,1)
 790         aad(i,2)=scalscp*aad(i,2)
 791         bad(i,1)=scalscp*bad(i,1)
 792         bad(i,2)=scalscp*bad(i,2)
 793       enddo
 794       wumb=1.0d0
 795       call rescale_weights(t_bath)
 796       if(me.eq.king.or..not.out1file)
 797      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 798      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 799      &  wturn4,wturn6
 800    22 format (/'Energy-term weights (scaled):'//
 801      & 'WSCC=   ',f10.6,' (SC-SC)'/
 802      & 'WSCP=   ',f10.6,' (SC-p)'/
 803      & 'WELEC=  ',f10.6,' (p-p electr)'/
 804      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 805      & 'WBOND=  ',f10.6,' (stretching)'/
 806      & 'WANG=   ',f10.6,' (bending)'/
 807      & 'WSCLOC= ',f10.6,' (SC local)'/
 808      & 'WTOR=   ',f10.6,' (torsional)'/
 809      & 'WTORD=  ',f10.6,' (double torsional)'/
 810      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 811      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 812      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 813      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 814      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 815      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 816      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 817      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 818      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 819       if(me.eq.king.or..not.out1file)
 820      & write (iout,*) "Reference temperature for weights calculation:",
 821      &  temp0
 822       call reada(weightcard,"D0CM",d0cm,3.78d0)
 823       call reada(weightcard,"AKCM",akcm,15.1d0)
 824       call reada(weightcard,"AKTH",akth,11.0d0)
 825       call reada(weightcard,"AKCT",akct,12.0d0)
 826       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 827       call reada(weightcard,"V2SS",v2ss,7.61d0)
 828       call reada(weightcard,"V3SS",v3ss,13.7d0)
 829       call reada(weightcard,"EBR",ebr,-5.50D0)
 830       call reada(weightcard,"ATRISS",atriss,0.301D0)
 831       call reada(weightcard,"BTRISS",btriss,0.021D0)
 832       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 833       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 834       write (iout,*) "ATRISS=", atriss
 835       write (iout,*) "BTRISS=", btriss
 836       write (iout,*) "CTRISS=", ctriss
 837       write (iout,*) "DTRISS=", dtriss
 838       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 839       do i=1,maxres
 840         dyn_ss_mask(i)=.false.
 841       enddo
 842       do i=1,maxres-1
 843         do j=i+1,maxres
 844           dyn_ssbond_ij(i,j)=1.0d300
 845         enddo
 846       enddo
 847       call reada(weightcard,"HT",Ht,0.0D0)
 848       if (dyn_ss) then
 849         ss_depth=ebr/wsc-0.25*eps(1,1)
 850         Ht=Ht/wsc-0.25*eps(1,1)
 851         akcm=akcm*wstrain/wsc
 852         akth=akth*wstrain/wsc
 853         akct=akct*wstrain/wsc
 854         v1ss=v1ss*wstrain/wsc
 855         v2ss=v2ss*wstrain/wsc
 856         v3ss=v3ss*wstrain/wsc
 857       else
 858         if (wstrain.ne.0.0) then
 859          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 860         else
 861           ss_depth=0.0
 862         endif
 863       endif
 864       write (iout,*) "wshield,", wshield
 865       if(me.eq.king.or..not.out1file) then
 866        write (iout,*) "Parameters of the SS-bond potential:"
 867        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 868      & " AKCT",akct
 869        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 870        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 871        write (iout,*)" HT",Ht
 872        print *,'indpdb=',indpdb,' pdbref=',pdbref
 873       endif
 874       if (indpdb.gt.0 .or. pdbref) then
 875         read(inp,'(a)') pdbfile
 876         if(me.eq.king.or..not.out1file)
 877      &   write (iout,'(2a)') 'PDB data will be read from file ',
 878      &   pdbfile(:ilen(pdbfile))
 879         open(ipdbin,file=pdbfile,status='old',err=33)
 880         goto 34
 881   33    write (iout,'(a)') 'Error opening PDB file.'
 882         stop
 883   34    continue
 884 c        write (iout,*) 'Begin reading pdb data'
 885 c        call flush(iout)
 886         call readpdb
 887 c        write (iout,*) 'Finished reading pdb data'
 888 c        call flush(iout)
 889         if(me.eq.king.or..not.out1file)
 890      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 891      &   ' nstart_sup=',nstart_sup
 892         do i=1,nres
 893           itype_pdb(i)=itype(i)
 894         enddo
 895         close (ipdbin)
 896         nnt=nstart_sup
 897         nct=nstart_sup+nsup-1
 898         call contact(.false.,ncont_ref,icont_ref,co)
 899
 900         if (sideadd) then
 901          if(me.eq.king.or..not.out1file)
 902      &    write(iout,*)'Adding sidechains'
 903          maxsi=1000
 904          do i=2,nres-1
 905           iti=itype(i)
 906           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 907             nsi=0
 908             fail=.true.
 909             do while (fail.and.nsi.le.maxsi)
 910               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 911               nsi=nsi+1
 912             enddo
 913             if(fail) write(iout,*)'Adding sidechain failed for res ',
 914      &              i,' after ',nsi,' trials'
 915           endif
 916          enddo
 917         endif
 918       endif
 919       if (indpdb.eq.0) then
 920 C Read sequence if not taken from the pdb file.
 921         read (inp,*) nres
 922 c        print *,'nres=',nres
 923         if (iscode.gt.0) then
 924           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 925         else
 926           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 927         endif
 928 C Convert sequence to numeric code
 929         do i=1,nres
 930           itype(i)=rescode(i,sequence(i),iscode)
 931         enddo
 932 C Assign initial virtual bond lengths
 933         do i=2,nres
 934           vbld(i)=vbl
 935           vbld_inv(i)=vblinv
 936         enddo
 937         do i=2,nres-1
 938           vbld(i+nres)=dsc(iabs(itype(i)))
 939           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 940 c          write (iout,*) "i",i," itype",itype(i),
 941 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 942         enddo
 943       endif
 944 c      print *,nres
 945 c      print '(20i4)',(itype(i),i=1,nres)
 946       do i=1,nres
 947 #ifdef PROCOR
 948         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 949 #else
 950         if (itype(i).eq.ntyp1) then
 951 #endif
 952           itel(i)=0
 953 #ifdef PROCOR
 954         else if (iabs(itype(i+1)).ne.20) then
 955 #else
 956         else if (iabs(itype(i)).ne.20) then
 957 #endif
 958           itel(i)=1
 959         else
 960           itel(i)=2
 961         endif
 962       enddo
 963       if(me.eq.king.or..not.out1file)then
 964        write (iout,*) "ITEL"
 965        do i=1,nres-1
 966          write (iout,*) i,itype(i),itel(i)
 967        enddo
 968        print *,'Call Read_Bridge.'
 969       endif
 970       nnt=1
 971       nct=nres
 972 cd      print *,'NNT=',NNT,' NCT=',NCT
 973       if (itype(1).eq.ntyp1) nnt=2
 974       if (itype(nres).eq.ntyp1) nct=nct-1
 975       if (pdbref) then
 976         if(me.eq.king.or..not.out1file)
 977      &   write (iout,'(a,i3)') 'nsup=',nsup
 978         nstart_seq=nnt
 979         if (nsup.le.(nct-nnt+1)) then
 980           do i=0,nct-nnt+1-nsup
 981             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 982               nstart_seq=nnt+i
 983               goto 111
 984             endif
 985           enddo
 986           write (iout,'(a)')
 987      &            'Error - sequences to be superposed do not match.'
 988           stop
 989         else
 990           do i=0,nsup-(nct-nnt+1)
 991             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 992      &      then
 993               nstart_sup=nstart_sup+i
 994               nsup=nct-nnt+1
 995               goto 111
 996             endif
 997           enddo
 998           write (iout,'(a)')
 999      &            'Error - sequences to be superposed do not match.'
1000         endif
1001   111   continue
1002         if (nsup.eq.0) nsup=nct-nnt
1003         if (nstart_sup.eq.0) nstart_sup=nnt
1004         if (nstart_seq.eq.0) nstart_seq=nnt
1005         if(me.eq.king.or..not.out1file)
1006      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1007      &                 ' nstart_seq=',nstart_seq
1008       endif
1009       call read_bridge
1010 C 8/13/98 Set limits to generating the dihedral angles
1011       do i=1,nres
1012         phibound(1,i)=-pi
1013         phibound(2,i)=pi
1014       enddo
1015       read (inp,*) ndih_constr
1016       write (iout,*) "ndih_constr",ndih_constr
1017       if (ndih_constr.gt.0) then
1018         raw_psipred=.false.
1019 C        read (inp,*) ftors
1020         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
1021      &  i=1,ndih_constr)
1022 #ifdef MPI
1023         if(me.eq.king.or..not.out1file)then
1024 #endif
1025          write (iout,*)
1026      &   'There are',ndih_constr,' restraints on gamma angles.'
1027          do i=1,ndih_constr
1028           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
1029      &    ftors(i)
1030          enddo
1031 #ifdef MPI
1032         endif
1033 #endif
1034         do i=1,ndih_constr
1035           phi0(i)=deg2rad*phi0(i)
1036           drange(i)=deg2rad*drange(i)
1037         enddo
1038 C        if(me.eq.king.or..not.out1file)
1039 C     &   write (iout,*) 'FTORS',ftors
1040         do i=1,ndih_constr
1041           ii = idih_constr(i)
1042           phibound(1,ii) = phi0(i)-drange(i)
1043           phibound(2,ii) = phi0(i)+drange(i)
1044         enddo
1045 #ifdef MPI
1046         if (me.eq.king .or. .not.out1file) then
1047 #endif
1048          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
1049          do i=1,nres
1050            write (iout,'(a3,i5,2f10.1)')
1051      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1052          enddo
1053 #ifdef MP
1054         endif
1055 #endif
1056       else if (ndih_constr.lt.0) then
1057         raw_psipred=.true.
1058         call card_concat(weightcard)
1059         call reada(weightcard,"PHIHEL",phihel,50.0D0)
1060         call reada(weightcard,"PHIBET",phibet,180.0D0)
1061         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1062         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1063         call reada(weightcard,"WDIHC",wdihc,0.591D0)
1064         write (iout,*) "Weight of dihedral angle restraints",wdihc
1065         read(inp,'(9x,3f7.3)')
1066      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1067         write (iout,*) "The secprob array"
1068         do i=nnt,nct
1069           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1070         enddo
1071         ndih_constr=0
1072         do i=nnt+3,nct
1073           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
1074      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
1075             ndih_constr=ndih_constr+1
1076             idih_constr(ndih_constr)=i
1077             sumv=0.0d0
1078             do j=1,3
1079               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1080               sumv=sumv+vpsipred(j,ndih_constr)
1081             enddo
1082             do j=1,3
1083               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1084             enddo
1085             phibound(1,ndih_constr)=phihel*deg2rad
1086             phibound(2,ndih_constr)=phibet*deg2rad
1087             sdihed(1,ndih_constr)=sigmahel*deg2rad
1088             sdihed(2,ndih_constr)=sigmabet*deg2rad
1089           endif
1090         enddo
1091 #ifdef MPI
1092         if(me.eq.king.or..not.out1file)then
1093 #endif
1094          write (iout,*)
1095      &   'There are',ndih_constr,
1096      &   ' bimodal restraints on gamma angles.'
1097          do i=1,ndih_constr
1098            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1099      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1100      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1101      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1102      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1103      &      (vpsipred(j,i),j=1,3)
1104          enddo
1105 #ifdef MPI
1106         endif
1107 #endif
1108 c Raw psipred input
1109       endif
1110 C first setting the theta boundaries to 0 to pi
1111 C this mean that there is no energy penalty for any angle occuring this can be applied
1112 C for generate random conformation but is not implemented in this way
1113 C      do i=1,nres
1114 C        thetabound(1,i)=0
1115 C        thetabound(2,i)=pi
1116 C      enddo
1117 C begin reading theta constrains this is quartic constrains allowing to
1118 C have smooth second derivative
1119       if (with_theta_constr) then
1120 C with_theta_constr is keyword allowing for occurance of theta constrains
1121       read (inp,*) ntheta_constr
1122 C ntheta_constr is the number of theta constrains
1123       if (ntheta_constr.gt.0) then
1124 C        read (inp,*) ftors
1125         read (inp,*) (itheta_constr(i),theta_constr0(i),
1126      &  theta_drange(i),for_thet_constr(i),
1127      &  i=1,ntheta_constr)
1128 C the above code reads from 1 to ntheta_constr
1129 C itheta_constr(i) residue i for which is theta_constr
1130 C theta_constr0 the global minimum value
1131 C theta_drange is range for which there is no energy penalty
1132 C for_thet_constr is the force constant for quartic energy penalty
1133 C E=k*x**4
1134         if(me.eq.king.or..not.out1file)then
1135          write (iout,*)
1136      &   'There are',ntheta_constr,' constraints on phi angles.'
1137          do i=1,ntheta_constr
1138           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1139      &    theta_drange(i),
1140      &    for_thet_constr(i)
1141          enddo
1142         endif
1143         do i=1,ntheta_constr
1144           theta_constr0(i)=deg2rad*theta_constr0(i)
1145           theta_drange(i)=deg2rad*theta_drange(i)
1146         enddo
1147 C        if(me.eq.king.or..not.out1file)
1148 C     &   write (iout,*) 'FTORS',ftors
1149 C        do i=1,ntheta_constr
1150 C          ii = itheta_constr(i)
1151 C          thetabound(1,ii) = phi0(i)-drange(i)
1152 C          thetabound(2,ii) = phi0(i)+drange(i)
1153 C        enddo
1154       endif ! ntheta_constr.gt.0
1155       endif! with_theta_constr
1156 C
1157 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1158 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1159 c--- Zscore rms -------
1160       if (nz_start.eq.0) nz_start=nnt
1161       if (nz_end.eq.0 .and. nsup.gt.0) then
1162         nz_end=nnt+nsup-1
1163       else if (nz_end.eq.0) then
1164         nz_end=nct
1165       endif
1166       if(me.eq.king.or..not.out1file)then
1167        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1168        write (iout,*) 'IZ_SC=',iz_sc
1169       endif
1170 c----------------------
1171       call init_int_table
1172       if (refstr) then
1173         if (.not.pdbref) then
1174           call read_angles(inp,*38)
1175           goto 39
1176    38     write (iout,'(a)') 'Error reading reference structure.'
1177 #ifdef MPI
1178           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1179           stop 'Error reading reference structure'
1180 #endif
1181    39     call chainbuild_extconf
1182           call setup_var
1183 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1184           nstart_sup=nnt
1185           nstart_seq=nnt
1186           nsup=nct-nnt+1
1187           kkk=1
1188           do i=1,2*nres
1189             do j=1,3
1190               cref(j,i,kkk)=c(j,i)
1191             enddo
1192           enddo
1193           call contact(.true.,ncont_ref,icont_ref,co)
1194         endif
1195         endif
1196 c        print *, "A TU"
1197         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1198         call flush(iout)
1199         if (constr_dist.gt.0) call read_dist_constr
1200         write (iout,*) "After read_dist_constr nhpb",nhpb
1201         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1202         call hpb_partition
1203         call NMRpeak_partition
1204         if(me.eq.king.or..not.out1file)
1205      &   write (iout,*) 'Contact order:',co
1206         if (pdbref) then
1207         if(me.eq.king.or..not.out1file)
1208      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1209         do i=1,ncont_ref
1210           do j=1,2
1211             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1212           enddo
1213           if(me.eq.king.or..not.out1file)
1214      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1215      &     icont_ref(1,i),' ',
1216      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1217         enddo
1218         endif
1219       write (iout,*) "calling read_saxs_consrtr",nsaxs
1220       if (nsaxs.gt.0) call read_saxs_constr
1221
1222 C      endif
1223       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1224      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1225      &    modecalc.ne.10) then
1226 C If input structure hasn't been supplied from the PDB file read or generate
1227 C initial geometry.
1228         if (iranconf.eq.0 .and. .not. extconf) then
1229           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1230      &     write (iout,'(a)') 'Initial geometry will be read in.'
1231           if (read_cart) then
1232             read(inp,'(8f10.5)',end=36,err=36)
1233      &       ((c(l,k),l=1,3),k=1,nres),
1234      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1235             write (iout,*) "Exit READ_CART"
1236 c            write (iout,'(8f10.5)')
1237 c     &       ((c(l,k),l=1,3),k=1,nres),
1238 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1239             call cartprint
1240             do i=1,nres-1
1241               do j=1,3
1242                 dc(j,i)=c(j,i+1)-c(j,i)
1243 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1244               enddo
1245             enddo
1246             do i=nnt,nct
1247               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1248                 do j=1,3
1249                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1250 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1251                 enddo
1252               else
1253                 do j=1,3
1254                   dc(j,i+nres)=0.0d0
1255 c                  dc_norm(j,i+nres)=0.0d0
1256                 enddo
1257               endif
1258             enddo
1259             call int_from_cart1(.true.)
1260             write (iout,*) "Finish INT_TO_CART"
1261 c      write (iout,*) "DC"
1262 c      do i=1,nres
1263 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1264 c     &    (dc(j,i+nres),j=1,3)
1265 c      enddo
1266 c            call cartprint
1267 c            call setup_var
1268 c            return
1269           else
1270             call read_angles(inp,*36)
1271             call chainbuild_extconf
1272           endif
1273           goto 37
1274    36     write (iout,'(a)') 'Error reading angle file.'
1275 #ifdef MPI
1276           call mpi_finalize( MPI_COMM_WORLD,IERR )
1277 #endif
1278           stop 'Error reading angle file.'
1279    37     continue
1280         else if (extconf) then
1281          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1282      &    write (iout,'(a)') 'Extended chain initial geometry.'
1283          do i=3,nres
1284           theta(i)=90d0*deg2rad
1285          enddo
1286          do i=4,nres
1287           phi(i)=180d0*deg2rad
1288          enddo
1289          do i=2,nres-1
1290           alph(i)=110d0*deg2rad
1291          enddo
1292          do i=2,nres-1
1293           omeg(i)=-120d0*deg2rad
1294           if (itype(i).le.0) omeg(i)=-omeg(i)
1295          enddo
1296         call chainbuild_extconf
1297         else
1298           if(me.eq.king.or..not.out1file)
1299      &     write (iout,'(a)') 'Random-generated initial geometry.'
1300
1301
1302 #ifdef MPI
1303           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1304      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1305 #endif
1306             do itrial=1,100
1307               itmp=1
1308               call gen_rand_conf(itmp,*30)
1309               goto 40
1310    30         write (iout,*) 'Failed to generate random conformation',
1311      &          ', itrial=',itrial
1312               write (*,*) 'Processor:',me,
1313      &          ' Failed to generate random conformation',
1314      &          ' itrial=',itrial
1315               call intout
1316
1317 #ifdef AIX
1318               call flush_(iout)
1319 #else
1320               call flush(iout)
1321 #endif
1322             enddo
1323             write (iout,'(a,i3,a)') 'Processor:',me,
1324      &        ' error in generating random conformation.'
1325             write (*,'(a,i3,a)') 'Processor:',me,
1326      &        ' error in generating random conformation.'
1327             call flush(iout)
1328 #ifdef MPI
1329             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1330    40       continue
1331           endif
1332 #else
1333           write (*,'(a)')
1334      &      ' error in generating random conformation.'
1335           stop
1336    40     continue
1337 #endif
1338         endif
1339       elseif (modecalc.eq.4) then
1340         read (inp,'(a)') intinname
1341         open (intin,file=intinname,status='old',err=333)
1342         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1343      &  write (iout,'(a)') 'intinname',intinname
1344         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1345         goto 334
1346   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1347 #ifdef MPI
1348         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1349 #endif
1350         stop 'Error opening angle file.'
1351   334   continue
1352
1353       endif
1354 C Generate distance constraints, if the PDB structure is to be regularized.
1355       if (nthread.gt.0) then
1356         call read_threadbase
1357       endif
1358       call setup_var
1359       if (me.eq.king .or. .not. out1file)
1360      & call intout
1361       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1362         write (iout,'(/a,i3,a)')
1363      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1364         write (iout,'(20i4)') (iss(i),i=1,ns)
1365        if (dyn_ss) then
1366           write(iout,*)"Running with dynamic disulfide-bond formation"
1367        else
1368         write (iout,'(/a/)') 'Pre-formed links are:'
1369         do i=1,nss
1370           i1=ihpb(i)-nres
1371           i2=jhpb(i)-nres
1372           it1=itype(i1)
1373           it2=itype(i2)
1374           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1375      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1376      &    ebr,forcon(i)
1377         enddo
1378         write (iout,'(a)')
1379        endif
1380       endif
1381       if (ns.gt.0.and.dyn_ss) then
1382           do i=nss+1,nhpb
1383             ihpb(i-nss)=ihpb(i)
1384             jhpb(i-nss)=jhpb(i)
1385             forcon(i-nss)=forcon(i)
1386             dhpb(i-nss)=dhpb(i)
1387           enddo
1388           nhpb=nhpb-nss
1389           nss=0
1390           call hpb_partition
1391           do i=1,ns
1392             dyn_ss_mask(iss(i))=.true.
1393           enddo
1394       endif
1395 c --- test
1396 c      call cartprint
1397 c      write (iout,*) "DC"
1398 c      do i=1,nres
1399 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1400 c     &    (dc(j,i+nres),j=1,3)
1401 c      enddo
1402       if (i2ndstr.gt.0) call secstrp2dihc
1403 c      call geom_to_var(nvar,x)
1404 c      call etotal(energia(0))
1405 c      call enerprint(energia(0))
1406 c      call briefout(0,etot)
1407 c      stop
1408 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1409 cd    write (iout,'(a)') 'Variable list:'
1410 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1411 #ifdef MPI
1412       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1413      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1414      &  'Processor',myrank,': end reading molecular data.'
1415 #endif
1416 c      print *,"A TU?"
1417       return
1418       end
1419 c--------------------------------------------------------------------------
1420       logical function seq_comp(itypea,itypeb,length)
1421       implicit none
1422       integer length,itypea(length),itypeb(length)
1423       integer i
1424       do i=1,length
1425         if (itypea(i).ne.itypeb(i)) then
1426           seq_comp=.false.
1427           return
1428         endif
1429       enddo
1430       seq_comp=.true.
1431       return
1432       end
1433 c-----------------------------------------------------------------------------
1434       subroutine read_bridge
1435 C Read information about disulfide bridges.
1436       implicit real*8 (a-h,o-z)
1437       include 'DIMENSIONS'
1438 #ifdef MPI
1439       include 'mpif.h'
1440 #endif
1441       include 'COMMON.IOUNITS'
1442       include 'COMMON.GEO'
1443       include 'COMMON.VAR'
1444       include 'COMMON.INTERACT'
1445       include 'COMMON.LOCAL'
1446       include 'COMMON.NAMES'
1447       include 'COMMON.CHAIN'
1448       include 'COMMON.FFIELD'
1449       include 'COMMON.SBRIDGE'
1450       include 'COMMON.HEADER'
1451       include 'COMMON.CONTROL'
1452       include 'COMMON.DBASE'
1453       include 'COMMON.THREAD'
1454       include 'COMMON.TIME1'
1455       include 'COMMON.SETUP'
1456 C Read bridging residues.
1457       read (inp,*) ns,(iss(i),i=1,ns)
1458       print *,'ns=',ns
1459       if(me.eq.king.or..not.out1file)
1460      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1461 C Check whether the specified bridging residues are cystines.
1462       do i=1,ns
1463         if (itype(iss(i)).ne.1) then
1464           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1465      &   'Do you REALLY think that the residue ',
1466      &    restyp(itype(iss(i))),i,
1467      &   ' can form a disulfide bridge?!!!'
1468           write (*,'(2a,i3,a)')
1469      &   'Do you REALLY think that the residue ',
1470      &    restyp(itype(iss(i))),i,
1471      &   ' can form a disulfide bridge?!!!'
1472 #ifdef MPI
1473          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1474          stop
1475 #endif
1476         endif
1477       enddo
1478 C Read preformed bridges.
1479       if (ns.gt.0) then
1480       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1481       if(fg_rank.eq.0)
1482      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1483       if (nss.gt.0) then
1484         nhpb=nss
1485 C Check if the residues involved in bridges are in the specified list of
1486 C bridging residues.
1487         do i=1,nss
1488           do j=1,i-1
1489             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1490      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1491               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1492      &      ' contains residues present in other pairs.'
1493               write (*,'(a,i3,a)') 'Disulfide pair',i,
1494      &      ' contains residues present in other pairs.'
1495 #ifdef MPI
1496               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1497               stop
1498 #endif
1499             endif
1500           enddo
1501           do j=1,ns
1502             if (ihpb(i).eq.iss(j)) goto 10
1503           enddo
1504           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1505    10     continue
1506           do j=1,ns
1507             if (jhpb(i).eq.iss(j)) goto 20
1508           enddo
1509           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1510    20     continue
1511           dhpb(i)=dbr
1512           forcon(i)=fbr
1513         enddo
1514         do i=1,nss
1515           ihpb(i)=ihpb(i)+nres
1516           jhpb(i)=jhpb(i)+nres
1517         enddo
1518       endif
1519       endif
1520       return
1521       end
1522 c----------------------------------------------------------------------------
1523       subroutine read_x(kanal,*)
1524       implicit real*8 (a-h,o-z)
1525       include 'DIMENSIONS'
1526       include 'COMMON.GEO'
1527       include 'COMMON.VAR'
1528       include 'COMMON.CHAIN'
1529       include 'COMMON.IOUNITS'
1530       include 'COMMON.CONTROL'
1531       include 'COMMON.LOCAL'
1532       include 'COMMON.INTERACT'
1533 c Read coordinates from input
1534 c
1535       read(kanal,'(8f10.5)',end=10,err=10)
1536      &  ((c(l,k),l=1,3),k=1,nres),
1537      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1538       do j=1,3
1539         c(j,nres+1)=c(j,1)
1540         c(j,2*nres)=c(j,nres)
1541       enddo
1542       call int_from_cart1(.false.)
1543       do i=1,nres-1
1544         do j=1,3
1545           dc(j,i)=c(j,i+1)-c(j,i)
1546           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1547         enddo
1548       enddo
1549       do i=nnt,nct
1550         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1551           do j=1,3
1552             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1553             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1554           enddo
1555         endif
1556       enddo
1557
1558       return
1559    10 return1
1560       end
1561 c----------------------------------------------------------------------------
1562       subroutine read_threadbase
1563       implicit real*8 (a-h,o-z)
1564       include 'DIMENSIONS'
1565       include 'COMMON.IOUNITS'
1566       include 'COMMON.GEO'
1567       include 'COMMON.VAR'
1568       include 'COMMON.INTERACT'
1569       include 'COMMON.LOCAL'
1570       include 'COMMON.NAMES'
1571       include 'COMMON.CHAIN'
1572       include 'COMMON.FFIELD'
1573       include 'COMMON.SBRIDGE'
1574       include 'COMMON.HEADER'
1575       include 'COMMON.CONTROL'
1576       include 'COMMON.DBASE'
1577       include 'COMMON.THREAD'
1578       include 'COMMON.TIME1'
1579 C Read pattern database for threading.
1580       read (icbase,*) nseq
1581       do i=1,nseq
1582         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1583      &   nres_base(2,i),nres_base(3,i)
1584         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1585      &   nres_base(1,i))
1586 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1587 c    &   nres_base(2,i),nres_base(3,i)
1588 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1589 c    &   nres_base(1,i))
1590       enddo
1591       close (icbase)
1592       if (weidis.eq.0.0D0) weidis=0.1D0
1593       do i=nnt,nct
1594         do j=i+2,nct
1595           nhpb=nhpb+1
1596           ihpb(nhpb)=i
1597           jhpb(nhpb)=j
1598           forcon(nhpb)=weidis
1599         enddo
1600       enddo
1601       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1602       write (iout,'(a,i5)') 'nexcl: ',nexcl
1603       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1604       return
1605       end
1606 c------------------------------------------------------------------------------
1607       subroutine setup_var
1608       implicit real*8 (a-h,o-z)
1609       include 'DIMENSIONS'
1610       include 'COMMON.IOUNITS'
1611       include 'COMMON.GEO'
1612       include 'COMMON.VAR'
1613       include 'COMMON.INTERACT'
1614       include 'COMMON.LOCAL'
1615       include 'COMMON.NAMES'
1616       include 'COMMON.CHAIN'
1617       include 'COMMON.FFIELD'
1618       include 'COMMON.SBRIDGE'
1619       include 'COMMON.HEADER'
1620       include 'COMMON.CONTROL'
1621       include 'COMMON.DBASE'
1622       include 'COMMON.THREAD'
1623       include 'COMMON.TIME1'
1624 C Set up variable list.
1625       ntheta=nres-2
1626       nphi=nres-3
1627       nvar=ntheta+nphi
1628       nside=0
1629       write (iout,*) "SETUP_VAR ialph"
1630       do i=2,nres-1
1631         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1632           nside=nside+1
1633           ialph(i,1)=nvar+nside
1634           ialph(nside,2)=i
1635         endif
1636       enddo
1637       if (indphi.gt.0) then
1638         nvar=nphi
1639       else if (indback.gt.0) then
1640         nvar=nphi+ntheta
1641       else
1642         nvar=nvar+2*nside
1643       endif
1644       write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1645       return
1646       end
1647 c----------------------------------------------------------------------------
1648       subroutine gen_dist_constr
1649 C Generate CA distance constraints.
1650       implicit real*8 (a-h,o-z)
1651       include 'DIMENSIONS'
1652       include 'COMMON.IOUNITS'
1653       include 'COMMON.GEO'
1654       include 'COMMON.VAR'
1655       include 'COMMON.INTERACT'
1656       include 'COMMON.LOCAL'
1657       include 'COMMON.NAMES'
1658       include 'COMMON.CHAIN'
1659       include 'COMMON.FFIELD'
1660       include 'COMMON.SBRIDGE'
1661       include 'COMMON.HEADER'
1662       include 'COMMON.CONTROL'
1663       include 'COMMON.DBASE'
1664       include 'COMMON.THREAD'
1665       include 'COMMON.TIME1'
1666       dimension itype_pdb(maxres)
1667       common /pizda/ itype_pdb
1668       character*2 iden
1669 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1670 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1671 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1672 cd     & ' nsup',nsup
1673       do i=nstart_sup,nstart_sup+nsup-1
1674 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1675 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1676         do j=i+2,nstart_sup+nsup-1
1677           nhpb=nhpb+1
1678           ihpb(nhpb)=i+nstart_seq-nstart_sup
1679           jhpb(nhpb)=j+nstart_seq-nstart_sup
1680           forcon(nhpb)=weidis
1681           dhpb(nhpb)=dist(i,j)
1682         enddo
1683       enddo
1684 cd      write (iout,'(a)') 'Distance constraints:'
1685 cd      do i=nss+1,nhpb
1686 cd        ii=ihpb(i)
1687 cd        jj=jhpb(i)
1688 cd        iden='CA'
1689 cd        if (ii.gt.nres) then
1690 cd          iden='SC'
1691 cd          ii=ii-nres
1692 cd          jj=jj-nres
1693 cd        endif
1694 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1695 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1696 cd     &  dhpb(i),forcon(i)
1697 cd      enddo
1698       return
1699       end
1700 c----------------------------------------------------------------------------
1701       subroutine map_read
1702       implicit real*8 (a-h,o-z)
1703       include 'DIMENSIONS'
1704       include 'COMMON.MAP'
1705       include 'COMMON.IOUNITS'
1706       character*3 angid(4) /'THE','PHI','ALP','OME'/
1707       character*80 mapcard,ucase
1708       do imap=1,nmap
1709         read (inp,'(a)') mapcard
1710         mapcard=ucase(mapcard)
1711         if (index(mapcard,'PHI').gt.0) then
1712           kang(imap)=1
1713         else if (index(mapcard,'THE').gt.0) then
1714           kang(imap)=2
1715         else if (index(mapcard,'ALP').gt.0) then
1716           kang(imap)=3
1717         else if (index(mapcard,'OME').gt.0) then
1718           kang(imap)=4
1719         else
1720           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1721           stop 'Error - illegal variable spec in MAP card.'
1722         endif
1723         call readi (mapcard,'RES1',res1(imap),0)
1724         call readi (mapcard,'RES2',res2(imap),0)
1725         if (res1(imap).eq.0) then
1726           res1(imap)=res2(imap)
1727         else if (res2(imap).eq.0) then
1728           res2(imap)=res1(imap)
1729         endif
1730         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1731           write (iout,'(a)')
1732      &    'Error - illegal definition of variable group in MAP.'
1733           stop 'Error - illegal definition of variable group in MAP.'
1734         endif
1735         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1736         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1737         call readi(mapcard,'NSTEP',nstep(imap),0)
1738         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1739           write (iout,'(a)')
1740      &     'Illegal boundary and/or step size specification in MAP.'
1741           stop 'Illegal boundary and/or step size specification in MAP.'
1742         endif
1743       enddo ! imap
1744       return
1745       end
1746 c----------------------------------------------------------------------------
1747       subroutine csaread
1748       implicit real*8 (a-h,o-z)
1749       include 'DIMENSIONS'
1750       include 'COMMON.IOUNITS'
1751       include 'COMMON.GEO'
1752       include 'COMMON.CSA'
1753       include 'COMMON.BANK'
1754       include 'COMMON.CONTROL'
1755       character*80 ucase
1756       character*620 mcmcard
1757       call card_concat(mcmcard)
1758
1759       call readi(mcmcard,'NCONF',nconf,50)
1760       call readi(mcmcard,'NADD',nadd,0)
1761       call readi(mcmcard,'JSTART',jstart,1)
1762       call readi(mcmcard,'JEND',jend,1)
1763       call readi(mcmcard,'NSTMAX',nstmax,500000)
1764       call readi(mcmcard,'N0',n0,1)
1765       call readi(mcmcard,'N1',n1,6)
1766       call readi(mcmcard,'N2',n2,4)
1767       call readi(mcmcard,'N3',n3,0)
1768       call readi(mcmcard,'N4',n4,0)
1769       call readi(mcmcard,'N5',n5,0)
1770       call readi(mcmcard,'N6',n6,10)
1771       call readi(mcmcard,'N7',n7,0)
1772       call readi(mcmcard,'N8',n8,0)
1773       call readi(mcmcard,'N9',n9,0)
1774       call readi(mcmcard,'N14',n14,0)
1775       call readi(mcmcard,'N15',n15,0)
1776       call readi(mcmcard,'N16',n16,0)
1777       call readi(mcmcard,'N17',n17,0)
1778       call readi(mcmcard,'N18',n18,0)
1779
1780       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1781
1782       call readi(mcmcard,'NDIFF',ndiff,2)
1783       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1784       call readi(mcmcard,'IS1',is1,1)
1785       call readi(mcmcard,'IS2',is2,8)
1786       call readi(mcmcard,'NRAN0',nran0,4)
1787       call readi(mcmcard,'NRAN1',nran1,2)
1788       call readi(mcmcard,'IRR',irr,1)
1789       call readi(mcmcard,'NSEED',nseed,20)
1790       call readi(mcmcard,'NTOTAL',ntotal,10000)
1791       call reada(mcmcard,'CUT1',cut1,2.0d0)
1792       call reada(mcmcard,'CUT2',cut2,5.0d0)
1793       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1794       call readi(mcmcard,'ICMAX',icmax,3)
1795       call readi(mcmcard,'IRESTART',irestart,0)
1796 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1797       ntbankm=0
1798 c!bankt
1799       call reada(mcmcard,'DELE',dele,20.0d0)
1800       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1801       call readi(mcmcard,'IREF',iref,0)
1802       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1803       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1804       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1805       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1806       write (iout,*) "NCONF_IN",nconf_in
1807       return
1808       end
1809 c----------------------------------------------------------------------------
1810 cfmc      subroutine mcmfread
1811 cfmc      implicit real*8 (a-h,o-z)
1812 cfmc      include 'DIMENSIONS'
1813 cfmc      include 'COMMON.MCMF'
1814 cfmc      include 'COMMON.IOUNITS'
1815 cfmc      include 'COMMON.GEO'
1816 cfmc      character*80 ucase
1817 cfmc      character*620 mcmcard
1818 cfmc      call card_concat(mcmcard)
1819 cfmc
1820 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1821 cfmc      write(iout,*)'MAXRANT=',maxrant
1822 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1823 cfmc      write(iout,*)'MAXFAM=',maxfam
1824 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1825 cfmc      write(iout,*)'NNET1=',nnet1
1826 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1827 cfmc      write(iout,*)'NNET2=',nnet2
1828 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1829 cfmc      write(iout,*)'NNET3=',nnet3
1830 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1831 cfmc      write(iout,*)'ILASTT=',ilastt
1832 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1833 cfmc      write(iout,*)'MAXSTR=',maxstr
1834 cfmc      maxstr_f=maxstr/maxfam
1835 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1836 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1837 cfmc      write(iout,*)'NMCMF=',nmcmf
1838 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1839 cfmc      write(iout,*)'IFOCUS=',ifocus
1840 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1841 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1842 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1843 cfmc      write(iout,*)'INTPRT=',intprt
1844 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1845 cfmc      write(iout,*)'IPRT=',iprt
1846 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1847 cfmc      write(iout,*)'IMAXTR=',imaxtr
1848 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1849 cfmc      write(iout,*)'MAXEVEN=',maxeven
1850 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1851 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1852 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1853 cfmc      write(iout,*)'INIMIN=',inimin
1854 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1855 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1856 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1857 cfmc      write(iout,*)'NTHREAD=',nthread
1858 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1859 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1860 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1861 cfmc      write(iout,*)'MAXPERT=',maxpert
1862 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1863 cfmc      write(iout,*)'IRMSD=',irmsd
1864 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1865 cfmc      write(iout,*)'DENEMIN=',denemin
1866 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1867 cfmc      write(iout,*)'RCUT1S=',rcut1s
1868 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1869 cfmc      write(iout,*)'RCUT1E=',rcut1e
1870 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1871 cfmc      write(iout,*)'RCUT2S=',rcut2s
1872 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1873 cfmc      write(iout,*)'RCUT2E=',rcut2e
1874 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1875 cfmc      write(iout,*)'DPERT1=',d_pert1
1876 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1877 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1878 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1879 cfmc      write(iout,*)'DPERT2=',d_pert2
1880 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1881 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1882 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1883 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1884 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1885 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1886 cfmc      d_pert1=deg2rad*d_pert1
1887 cfmc      d_pert1a=deg2rad*d_pert1a
1888 cfmc      d_pert2=deg2rad*d_pert2
1889 cfmc      d_pert2a=deg2rad*d_pert2a
1890 cfmc      d_pert2b=deg2rad*d_pert2b
1891 cfmc      d_pert2c=deg2rad*d_pert2c
1892 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1893 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1894 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1895 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1896 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1897 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1898 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1899 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1900 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1901 cfmc      write(iout,*)'RCUTINI=',rcutini
1902 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1903 cfmc      write(iout,*)'GRAT=',grat
1904 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1905 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1906 cfmc
1907 cfmc      return
1908 cfmc      end
1909 c----------------------------------------------------------------------------
1910       subroutine mcmread
1911       implicit real*8 (a-h,o-z)
1912       include 'DIMENSIONS'
1913       include 'COMMON.MCM'
1914       include 'COMMON.MCE'
1915       include 'COMMON.IOUNITS'
1916       character*80 ucase
1917       character*320 mcmcard
1918       call card_concat(mcmcard)
1919       call readi(mcmcard,'MAXACC',maxacc,100)
1920       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1921       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1922       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1923       call readi(mcmcard,'MAXREPM',maxrepm,200)
1924       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1925       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1926       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1927       call reada(mcmcard,'E_UP',e_up,5.0D0)
1928       call reada(mcmcard,'DELTE',delte,0.1D0)
1929       call readi(mcmcard,'NSWEEP',nsweep,5)
1930       call readi(mcmcard,'NSTEPH',nsteph,0)
1931       call readi(mcmcard,'NSTEPC',nstepc,0)
1932       call reada(mcmcard,'TMIN',tmin,298.0D0)
1933       call reada(mcmcard,'TMAX',tmax,298.0D0)
1934       call readi(mcmcard,'NWINDOW',nwindow,0)
1935       call readi(mcmcard,'PRINT_MC',print_mc,0)
1936       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1937       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1938       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1939       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1940       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1941       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1942       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1943       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1944       if (nwindow.gt.0) then
1945         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1946         do i=1,nwindow
1947           winlen(i)=winend(i)-winstart(i)+1
1948         enddo
1949       endif
1950       if (tmax.lt.tmin) tmax=tmin
1951       if (tmax.eq.tmin) then
1952         nstepc=0
1953         nsteph=0
1954       endif
1955       if (nstepc.gt.0 .and. nsteph.gt.0) then
1956         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1957         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1958       endif
1959 C Probabilities of different move types
1960       sumpro_type(0)=0.0D0
1961       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1962       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1963       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1964       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1965       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1966       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1967       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1968       do i=1,MaxMoveType
1969         print *,'i',i,' sumprotype',sumpro_type(i)
1970         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1971         print *,'i',i,' sumprotype',sumpro_type(i)
1972       enddo
1973       return
1974       end
1975 c----------------------------------------------------------------------------
1976       subroutine read_minim
1977       implicit real*8 (a-h,o-z)
1978       include 'DIMENSIONS'
1979       include 'COMMON.MINIM'
1980       include 'COMMON.IOUNITS'
1981       character*80 ucase
1982       character*320 minimcard
1983       call card_concat(minimcard)
1984       call readi(minimcard,'MAXMIN',maxmin,2000)
1985       call readi(minimcard,'MAXFUN',maxfun,5000)
1986       call readi(minimcard,'MINMIN',minmin,maxmin)
1987       call readi(minimcard,'MINFUN',minfun,maxmin)
1988       call reada(minimcard,'TOLF',tolf,1.0D-2)
1989       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1990       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1991       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1992       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1993       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1994      &         'Options in energy minimization:'
1995       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1996      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1997      & 'MinMin:',MinMin,' MinFun:',MinFun,
1998      & ' TolF:',TolF,' RTolF:',RTolF
1999       return
2000       end
2001 c----------------------------------------------------------------------------
2002       subroutine read_angles(kanal,*)
2003       implicit real*8 (a-h,o-z)
2004       include 'DIMENSIONS'
2005       include 'COMMON.GEO'
2006       include 'COMMON.VAR'
2007       include 'COMMON.CHAIN'
2008       include 'COMMON.IOUNITS'
2009       include 'COMMON.CONTROL'
2010 c Read angles from input
2011 c
2012        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2013        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2014        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2015        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2016
2017        do i=1,nres
2018 c 9/7/01 avoid 180 deg valence angle
2019         if (theta(i).gt.179.99d0) theta(i)=179.99d0
2020 c
2021         theta(i)=deg2rad*theta(i)
2022         phi(i)=deg2rad*phi(i)
2023         alph(i)=deg2rad*alph(i)
2024         omeg(i)=deg2rad*omeg(i)
2025        enddo
2026       return
2027    10 return1
2028       end
2029 c----------------------------------------------------------------------------
2030       subroutine reada(rekord,lancuch,wartosc,default)
2031       implicit none
2032       character*(*) rekord,lancuch
2033       double precision wartosc,default
2034       integer ilen,iread
2035       external ilen
2036       iread=index(rekord,lancuch)
2037       if (iread.eq.0) then
2038         wartosc=default
2039         return
2040       endif
2041       iread=iread+ilen(lancuch)+1
2042       read (rekord(iread:),*,err=10,end=10) wartosc
2043       return
2044   10  wartosc=default
2045       return
2046       end
2047 c----------------------------------------------------------------------------
2048       subroutine readi(rekord,lancuch,wartosc,default)
2049       implicit none
2050       character*(*) rekord,lancuch
2051       integer wartosc,default
2052       integer ilen,iread
2053       external ilen
2054       iread=index(rekord,lancuch)
2055       if (iread.eq.0) then
2056         wartosc=default
2057         return
2058       endif
2059       iread=iread+ilen(lancuch)+1
2060       read (rekord(iread:),*,err=10,end=10) wartosc
2061       return
2062   10  wartosc=default
2063       return
2064       end
2065 c----------------------------------------------------------------------------
2066       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2067       implicit none
2068       integer dim,i
2069       integer tablica(dim),default
2070       character*(*) rekord,lancuch
2071       character*80 aux
2072       integer ilen,iread
2073       external ilen
2074       do i=1,dim
2075         tablica(i)=default
2076       enddo
2077       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2078       if (iread.eq.0) return
2079       iread=iread+ilen(lancuch)+1
2080       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2081    10 return
2082       end
2083 c----------------------------------------------------------------------------
2084       subroutine multreada(rekord,lancuch,tablica,dim,default)
2085       implicit none
2086       integer dim,i
2087       double precision tablica(dim),default
2088       character*(*) rekord,lancuch
2089       character*80 aux
2090       integer ilen,iread
2091       external ilen
2092       do i=1,dim
2093         tablica(i)=default
2094       enddo
2095       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2096       if (iread.eq.0) return
2097       iread=iread+ilen(lancuch)+1
2098       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2099    10 return
2100       end
2101 c----------------------------------------------------------------------------
2102       subroutine openunits
2103       implicit real*8 (a-h,o-z)
2104       include 'DIMENSIONS'
2105 #ifdef MPI
2106       include 'mpif.h'
2107       character*16 form,nodename
2108       integer nodelen
2109 #endif
2110       include 'COMMON.SETUP'
2111       include 'COMMON.IOUNITS'
2112       include 'COMMON.MD'
2113       include 'COMMON.CONTROL'
2114       integer lenpre,lenpot,ilen,lentmp
2115       external ilen
2116       character*3 out1file_text,ucase
2117       character*3 ll
2118       external ucase
2119       tmpdir=""
2120       out1file_text="YES"
2121 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2122       call getenv_loc("PREFIX",prefix)
2123       pref_orig = prefix
2124       call getenv_loc("POT",pot)
2125       call getenv_loc("DIRTMP",tmpdir)
2126       call getenv_loc("CURDIR",curdir)
2127       call getenv_loc("OUT1FILE",out1file_text)
2128 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2129       out1file_text=ucase(out1file_text)
2130       if (out1file_text(1:1).eq."Y") then
2131         out1file=.true.
2132       else
2133         out1file=fg_rank.gt.0
2134       endif
2135       lenpre=ilen(prefix)
2136       lenpot=ilen(pot)
2137       lentmp=ilen(tmpdir)
2138       if (lentmp.gt.0) then
2139           write (*,'(80(1h!))')
2140           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2141           write (*,'(80(1h!))')
2142           write (*,*)"All output files will be on node /tmp directory."
2143 #ifdef MPI
2144         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2145         if (me.eq.king) then
2146           write (*,*) "The master node is ",nodename
2147         else if (fg_rank.eq.0) then
2148           write (*,*) "I am the CG slave node ",nodename
2149         else
2150           write (*,*) "I am the FG slave node ",nodename
2151         endif
2152 #endif
2153         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2154         lenpre = lentmp+lenpre+1
2155       endif
2156       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2157 C Get the names and open the input files
2158 #if defined(WINIFL) || defined(WINPGI)
2159       open(1,file=pref_orig(:ilen(pref_orig))//
2160      &  '.inp',status='old',readonly,shared)
2161        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2162 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2163 C Get parameter filenames and open the parameter files.
2164       call getenv_loc('BONDPAR',bondname)
2165       open (ibond,file=bondname,status='old',readonly,shared)
2166       call getenv_loc('THETPAR',thetname)
2167       open (ithep,file=thetname,status='old',readonly,shared)
2168       call getenv_loc('ROTPAR',rotname)
2169       open (irotam,file=rotname,status='old',readonly,shared)
2170       call getenv_loc('TORPAR',torname)
2171       open (itorp,file=torname,status='old',readonly,shared)
2172       call getenv_loc('TORDPAR',tordname)
2173       open (itordp,file=tordname,status='old',readonly,shared)
2174       call getenv_loc('FOURIER',fouriername)
2175       open (ifourier,file=fouriername,status='old',readonly,shared)
2176       call getenv_loc('ELEPAR',elename)
2177       open (ielep,file=elename,status='old',readonly,shared)
2178       call getenv_loc('SIDEPAR',sidename)
2179       open (isidep,file=sidename,status='old',readonly,shared)
2180       call getenv_loc('LIPTRANPAR',liptranname)
2181       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2182 #elif (defined CRAY) || (defined AIX)
2183       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2184      &  action='read')
2185 c      print *,"Processor",myrank," opened file 1"
2186       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2187 c      print *,"Processor",myrank," opened file 9"
2188 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2189 C Get parameter filenames and open the parameter files.
2190       call getenv_loc('BONDPAR',bondname)
2191       open (ibond,file=bondname,status='old',action='read')
2192 c      print *,"Processor",myrank," opened file IBOND"
2193       call getenv_loc('THETPAR',thetname)
2194       open (ithep,file=thetname,status='old',action='read')
2195 #ifndef CRYST_THETA
2196       call getenv_loc('THETPARPDB',thetname_pdb)
2197       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2198 #endif
2199 c      print *,"Processor",myrank," opened file ITHEP"
2200       call getenv_loc('ROTPAR',rotname)
2201       open (irotam,file=rotname,status='old',action='read')
2202 #ifndef CRYST_SC
2203       call getenv_loc('ROTPARPDB',rotname_pdb)
2204       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2205 #endif
2206 c      print *,"Processor",myrank," opened file IROTAM"
2207       call getenv_loc('TORPAR',torname)
2208       open (itorp,file=torname,status='old',action='read')
2209 c      print *,"Processor",myrank," opened file ITORP"
2210       call getenv_loc('TORDPAR',tordname)
2211       open (itordp,file=tordname,status='old',action='read')
2212 c      print *,"Processor",myrank," opened file ITORDP"
2213       call getenv_loc('SCCORPAR',sccorname)
2214       open (isccor,file=sccorname,status='old',action='read')
2215 c      print *,"Processor",myrank," opened file ISCCOR"
2216       call getenv_loc('FOURIER',fouriername)
2217       open (ifourier,file=fouriername,status='old',action='read')
2218 c      print *,"Processor",myrank," opened file IFOURIER"
2219       call getenv_loc('ELEPAR',elename)
2220       open (ielep,file=elename,status='old',action='read')
2221 c      print *,"Processor",myrank," opened file IELEP"
2222       call getenv_loc('SIDEPAR',sidename)
2223       open (isidep,file=sidename,status='old',action='read')
2224       call getenv_loc('LIPTRANPAR',liptranname)
2225       open (iliptranpar,file=liptranname,status='old',action='read')
2226 c      print *,"Processor",myrank," opened file ISIDEP"
2227 c      print *,"Processor",myrank," opened parameter files"
2228 #elif (defined G77)
2229       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2230       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2231 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2232 C Get parameter filenames and open the parameter files.
2233       call getenv_loc('BONDPAR',bondname)
2234       open (ibond,file=bondname,status='old')
2235       call getenv_loc('THETPAR',thetname)
2236       open (ithep,file=thetname,status='old')
2237 #ifndef CRYST_THETA
2238       call getenv_loc('THETPARPDB',thetname_pdb)
2239       open (ithep_pdb,file=thetname_pdb,status='old')
2240 #endif
2241       call getenv_loc('ROTPAR',rotname)
2242       open (irotam,file=rotname,status='old')
2243 #ifndef CRYST_SC
2244       call getenv_loc('ROTPARPDB',rotname_pdb)
2245       open (irotam_pdb,file=rotname_pdb,status='old')
2246 #endif
2247       call getenv_loc('TORPAR',torname)
2248       open (itorp,file=torname,status='old')
2249       call getenv_loc('TORDPAR',tordname)
2250       open (itordp,file=tordname,status='old')
2251       call getenv_loc('SCCORPAR',sccorname)
2252       open (isccor,file=sccorname,status='old')
2253       call getenv_loc('FOURIER',fouriername)
2254       open (ifourier,file=fouriername,status='old')
2255       call getenv_loc('ELEPAR',elename)
2256       open (ielep,file=elename,status='old')
2257       call getenv_loc('SIDEPAR',sidename)
2258       open (isidep,file=sidename,status='old')
2259       call getenv_loc('LIPTRANPAR',liptranname)
2260       open (iliptranpar,file=liptranname,status='old')
2261 #else
2262       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2263      &  readonly)
2264        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2265 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2266 C Get parameter filenames and open the parameter files.
2267       call getenv_loc('BONDPAR',bondname)
2268       open (ibond,file=bondname,status='old',readonly)
2269       call getenv_loc('THETPAR',thetname)
2270       open (ithep,file=thetname,status='old',readonly)
2271       call getenv_loc('ROTPAR',rotname)
2272       open (irotam,file=rotname,status='old',readonly)
2273       call getenv_loc('TORPAR',torname)
2274       open (itorp,file=torname,status='old',readonly)
2275       call getenv_loc('TORDPAR',tordname)
2276       open (itordp,file=tordname,status='old',readonly)
2277       call getenv_loc('SCCORPAR',sccorname)
2278       open (isccor,file=sccorname,status='old',readonly)
2279 #ifndef CRYST_THETA
2280       call getenv_loc('THETPARPDB',thetname_pdb)
2281       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2282 #endif
2283       call getenv_loc('FOURIER',fouriername)
2284       open (ifourier,file=fouriername,status='old',readonly)
2285       call getenv_loc('ELEPAR',elename)
2286       open (ielep,file=elename,status='old',readonly)
2287       call getenv_loc('SIDEPAR',sidename)
2288       open (isidep,file=sidename,status='old',readonly)
2289       call getenv_loc('LIPTRANPAR',liptranname)
2290       open (iliptranpar,file=liptranname,status='old',action='read')
2291 #ifndef CRYST_SC
2292       call getenv_loc('ROTPARPDB',rotname_pdb)
2293       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2294 #endif
2295 #endif
2296 #ifdef TUBE
2297       call getenv_loc('TUBEPAR',tubename)
2298 #if defined(WINIFL) || defined(WINPGI)
2299       open (itube,file=tubename,status='old',readonly,shared)
2300 #elif (defined CRAY)  || (defined AIX)
2301       open (itube,file=tubename,status='old',action='read')
2302 #elif (defined G77)
2303       open (itube,file=tubename,status='old')
2304 #else
2305       open (itube,file=tubename,status='old',readonly)
2306 #endif
2307 #endif
2308 #ifndef OLDSCP
2309 C
2310 C 8/9/01 In the newest version SCp interaction constants are read from a file
2311 C Use -DOLDSCP to use hard-coded constants instead.
2312 C
2313       call getenv_loc('SCPPAR',scpname)
2314 #if defined(WINIFL) || defined(WINPGI)
2315       open (iscpp,file=scpname,status='old',readonly,shared)
2316 #elif (defined CRAY)  || (defined AIX)
2317       open (iscpp,file=scpname,status='old',action='read')
2318 #elif (defined G77)
2319       open (iscpp,file=scpname,status='old')
2320 #else
2321       open (iscpp,file=scpname,status='old',readonly)
2322 #endif
2323 #endif
2324       call getenv_loc('PATTERN',patname)
2325 #if defined(WINIFL) || defined(WINPGI)
2326       open (icbase,file=patname,status='old',readonly,shared)
2327 #elif (defined CRAY)  || (defined AIX)
2328       open (icbase,file=patname,status='old',action='read')
2329 #elif (defined G77)
2330       open (icbase,file=patname,status='old')
2331 #else
2332       open (icbase,file=patname,status='old',readonly)
2333 #endif
2334 #ifdef MPI
2335 C Open output file only for CG processes
2336 c      print *,"Processor",myrank," fg_rank",fg_rank
2337       if (fg_rank.eq.0) then
2338
2339       if (nodes.eq.1) then
2340         npos=3
2341       else
2342         npos = dlog10(dfloat(nodes-1))+1
2343       endif
2344       if (npos.lt.3) npos=3
2345       write (liczba,'(i1)') npos
2346       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2347      &  //')'
2348       write (liczba,form) me
2349       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2350      &  liczba(:ilen(liczba))
2351       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2352      &  //'.int'
2353       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2354      &  //'.pdb'
2355       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2356      &  liczba(:ilen(liczba))//'.mol2'
2357       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2358      &  liczba(:ilen(liczba))//'.stat'
2359       if (lentmp.gt.0)
2360      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2361      &      //liczba(:ilen(liczba))//'.stat')
2362       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2363      &  //'.rst'
2364       if(usampl) then
2365           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2366      & liczba(:ilen(liczba))//'.const'
2367       endif
2368
2369       endif
2370 #else
2371       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2372       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2373       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2374       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2375       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2376       if (lentmp.gt.0)
2377      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2378      &    //'.stat')
2379       rest2name=prefix(:ilen(prefix))//'.rst'
2380       if(usampl) then
2381          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2382       endif
2383 #endif
2384 #if defined(AIX) || defined(PGI) || defined(CRAY)
2385       if (me.eq.king .or. .not. out1file) then
2386          open(iout,file=outname,status='unknown')
2387       else
2388          open(iout,file="/dev/null",status="unknown")
2389       endif
2390 c#define DEBUG
2391 #ifdef DEBUG
2392       if (fg_rank.gt.0) then
2393         write (liczba,'(i3.3)') myrank/nfgtasks
2394         write (ll,'(bz,i3.3)') fg_rank
2395         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2396      &   status='unknown')
2397       endif
2398 #endif
2399 c#undef DEBUG
2400       if(me.eq.king) then
2401        open(igeom,file=intname,status='unknown',position='append')
2402        open(ipdb,file=pdbname,status='unknown')
2403        open(imol2,file=mol2name,status='unknown')
2404        open(istat,file=statname,status='unknown',position='append')
2405       else
2406 c1out       open(iout,file=outname,status='unknown')
2407       endif
2408 #else
2409       if (me.eq.king .or. .not.out1file)
2410      &    open(iout,file=outname,status='unknown')
2411 #ifdef DEBUG
2412       if (fg_rank.gt.0) then
2413         write (liczba,'(i3.3)') myrank/nfgtasks
2414         write (ll,'(bz,i3.3)') fg_rank
2415         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2416      &   status='unknown')
2417       endif
2418 #endif
2419       if(me.eq.king) then
2420        open(igeom,file=intname,status='unknown',access='append')
2421        open(ipdb,file=pdbname,status='unknown')
2422        open(imol2,file=mol2name,status='unknown')
2423        open(istat,file=statname,status='unknown',access='append')
2424       else
2425 c1out       open(iout,file=outname,status='unknown')
2426       endif
2427 #endif
2428       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2429       csa_seed=prefix(:lenpre)//'.CSA.seed'
2430       csa_history=prefix(:lenpre)//'.CSA.history'
2431       csa_bank=prefix(:lenpre)//'.CSA.bank'
2432       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2433       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2434       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2435 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2436       csa_int=prefix(:lenpre)//'.int'
2437       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2438       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2439       csa_in=prefix(:lenpre)//'.CSA.in'
2440 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2441 C Write file names
2442       if (me.eq.king)then
2443       write (iout,'(80(1h-))')
2444       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2445       write (iout,'(80(1h-))')
2446       write (iout,*) "Input file                      : ",
2447      &  pref_orig(:ilen(pref_orig))//'.inp'
2448       write (iout,*) "Output file                     : ",
2449      &  outname(:ilen(outname))
2450       write (iout,*)
2451       write (iout,*) "Sidechain potential file        : ",
2452      &  sidename(:ilen(sidename))
2453 #ifndef OLDSCP
2454       write (iout,*) "SCp potential file              : ",
2455      &  scpname(:ilen(scpname))
2456 #endif
2457       write (iout,*) "Electrostatic potential file    : ",
2458      &  elename(:ilen(elename))
2459       write (iout,*) "Cumulant coefficient file       : ",
2460      &  fouriername(:ilen(fouriername))
2461       write (iout,*) "Torsional parameter file        : ",
2462      &  torname(:ilen(torname))
2463       write (iout,*) "Double torsional parameter file : ",
2464      &  tordname(:ilen(tordname))
2465       write (iout,*) "SCCOR parameter file : ",
2466      &  sccorname(:ilen(sccorname))
2467       write (iout,*) "Bond & inertia constant file    : ",
2468      &  bondname(:ilen(bondname))
2469       write (iout,*) "Bending-torsion parameter file  : ",
2470      &  thetname(:ilen(thetname))
2471       write (iout,*) "Rotamer parameter file          : ",
2472      &  rotname(:ilen(rotname))
2473       write (iout,*) "Thetpdb parameter file          : ",
2474      &  thetname_pdb(:ilen(thetname_pdb))
2475       write (iout,*) "Threading database              : ",
2476      &  patname(:ilen(patname))
2477       if (lentmp.ne.0)
2478      &write (iout,*)" DIRTMP                          : ",
2479      &  tmpdir(:lentmp)
2480       write (iout,'(80(1h-))')
2481       endif
2482       return
2483       end
2484 c----------------------------------------------------------------------------
2485       subroutine card_concat(card)
2486       implicit real*8 (a-h,o-z)
2487       include 'DIMENSIONS'
2488       include 'COMMON.IOUNITS'
2489       character*(*) card
2490       character*80 karta,ucase
2491       external ilen
2492       read (inp,'(a)') karta
2493       karta=ucase(karta)
2494       card=' '
2495       do while (karta(80:80).eq.'&')
2496         card=card(:ilen(card)+1)//karta(:79)
2497         read (inp,'(a)') karta
2498         karta=ucase(karta)
2499       enddo
2500       card=card(:ilen(card)+1)//karta
2501       return
2502       end
2503 c----------------------------------------------------------------------------------
2504       subroutine readrst
2505       implicit real*8 (a-h,o-z)
2506       include 'DIMENSIONS'
2507       include 'COMMON.CHAIN'
2508       include 'COMMON.IOUNITS'
2509       include 'COMMON.MD'
2510       open(irest2,file=rest2name,status='unknown')
2511       read(irest2,*) totT,EK,potE,totE,t_bath
2512       totTafm=totT
2513       do i=1,2*nres
2514          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2515       enddo
2516       do i=1,2*nres
2517          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2518       enddo
2519       if(usampl) then
2520              read (irest2,*) iset
2521       endif
2522       close(irest2)
2523       return
2524       end
2525 c---------------------------------------------------------------------------------
2526       subroutine read_fragments
2527       implicit real*8 (a-h,o-z)
2528       include 'DIMENSIONS'
2529 #ifdef MPI
2530       include 'mpif.h'
2531 #endif
2532       include 'COMMON.SETUP'
2533       include 'COMMON.CHAIN'
2534       include 'COMMON.IOUNITS'
2535       include 'COMMON.MD'
2536       include 'COMMON.CONTROL'
2537       read(inp,*) nset,nfrag,npair,nfrag_back
2538       loc_qlike=(nfrag_back.lt.0)
2539       nfrag_back=iabs(nfrag_back)
2540       if(me.eq.king.or..not.out1file)
2541      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2542      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2543       do iset=1,nset
2544          read(inp,*) mset(iset)
2545        do i=1,nfrag
2546          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2547      &     qinfrag(i,iset)
2548          if(me.eq.king.or..not.out1file)
2549      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2550      &     ifrag(2,i,iset), qinfrag(i,iset)
2551        enddo
2552        do i=1,npair
2553         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2554      &    qinpair(i,iset)
2555         if(me.eq.king.or..not.out1file)
2556      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2557      &    ipair(2,i,iset), qinpair(i,iset)
2558        enddo
2559        if (loc_qlike) then
2560        do i=1,nfrag_back
2561         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2562      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2563      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2564      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2565         if(me.eq.king.or..not.out1file)
2566      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2567      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2568      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2569      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2570        enddo
2571        else
2572        do i=1,nfrag_back
2573         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2574      &     wfrag_back(3,i,iset),
2575      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2576         if(me.eq.king.or..not.out1file)
2577      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2578      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2579        enddo
2580        endif
2581       enddo
2582       return
2583       end
2584 C---------------------------------------------------------------------------
2585       subroutine read_afminp
2586             implicit real*8 (a-h,o-z)
2587       include 'DIMENSIONS'
2588 #ifdef MPI
2589       include 'mpif.h'
2590 #endif
2591       include 'COMMON.SETUP'
2592       include 'COMMON.CONTROL'
2593       include 'COMMON.CHAIN'
2594       include 'COMMON.IOUNITS'
2595       include 'COMMON.SBRIDGE'
2596       character*320 afmcard
2597       print *, "wchodze"
2598       call card_concat(afmcard)
2599       call readi(afmcard,"BEG",afmbeg,0)
2600       call readi(afmcard,"END",afmend,0)
2601       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2602       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2603       print *,'FORCE=' ,forceAFMconst
2604 CCCC NOW PROPERTIES FOR AFM
2605        distafminit=0.0d0
2606        do i=1,3
2607         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2608        enddo
2609         distafminit=dsqrt(distafminit)
2610         print *,'initdist',distafminit
2611       return
2612       end
2613 c-------------------------------------------------------------------------------
2614       subroutine read_saxs_constr
2615       implicit real*8 (a-h,o-z)
2616       include 'DIMENSIONS'
2617 #ifdef MPI
2618       include 'mpif.h'
2619 #endif
2620       include 'COMMON.SETUP'
2621       include 'COMMON.CONTROL'
2622       include 'COMMON.CHAIN'
2623       include 'COMMON.IOUNITS'
2624       include 'COMMON.SBRIDGE'
2625       double precision cm(3)
2626 c      read(inp,*) nsaxs
2627       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2628       call flush(iout)
2629       if (saxs_mode.eq.0) then
2630 c SAXS distance distribution
2631       do i=1,nsaxs
2632         read(inp,*) distsaxs(i),Psaxs(i)
2633       enddo
2634       Cnorm = 0.0d0
2635       do i=1,nsaxs
2636         Cnorm = Cnorm + Psaxs(i)
2637       enddo
2638       write (iout,*) "Cnorm",Cnorm
2639       do i=1,nsaxs
2640         Psaxs(i)=Psaxs(i)/Cnorm
2641       enddo
2642       write (iout,*) "Normalized distance distribution from SAXS"
2643       do i=1,nsaxs
2644         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2645       enddo
2646       Wsaxs0=0.0d0
2647       do i=1,nsaxs
2648         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2649       enddo
2650       write (iout,*) "Wsaxs0",Wsaxs0
2651       else
2652 c SAXS "spheres".
2653       do i=1,nsaxs
2654         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2655       enddo
2656       do j=1,3
2657         cm(j)=0.0d0
2658       enddo
2659       do i=1,nsaxs
2660         do j=1,3
2661           cm(j)=cm(j)+Csaxs(j,i)
2662         enddo
2663       enddo
2664       do j=1,3
2665         cm(j)=cm(j)/nsaxs
2666       enddo
2667       do i=1,nsaxs
2668         do j=1,3
2669           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2670         enddo
2671       enddo
2672       write (iout,*) "SAXS sphere coordinates"
2673       do i=1,nsaxs
2674         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2675       enddo
2676       endif
2677       return
2678       end
2679
2680 c-------------------------------------------------------------------------------
2681       subroutine read_dist_constr
2682       implicit real*8 (a-h,o-z)
2683       include 'DIMENSIONS'
2684 #ifdef MPI
2685       include 'mpif.h'
2686 #endif
2687       include 'COMMON.SETUP'
2688       include 'COMMON.CONTROL'
2689       include 'COMMON.CHAIN'
2690       include 'COMMON.IOUNITS'
2691       include 'COMMON.SBRIDGE'
2692       integer ifrag_(2,100),ipair_(2,1000)
2693       double precision wfrag_(100),wpair_(1000)
2694       character*5000 controlcard
2695       logical normalize,next
2696       integer restr_type
2697       double precision xlink(4,0:4) /
2698 c           a          b       c     sigma
2699      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2700      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2701      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2702      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2703      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2704 c      print *, "WCHODZE"
2705       write (iout,*) "Calling read_dist_constr"
2706 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2707 c      call flush(iout)
2708       next=.true.
2709
2710       npeak=0
2711       ipeak=0
2712       nhpb_peak=0
2713
2714       DO WHILE (next)
2715
2716       call card_concat(controlcard)
2717       next = index(controlcard,"NEXT").gt.0
2718       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2719       write (iout,*) "restr_type",restr_type
2720       call readi(controlcard,"NFRAG",nfrag_,0)
2721       call readi(controlcard,"NFRAG",nfrag_,0)
2722       call readi(controlcard,"NPAIR",npair_,0)
2723       call readi(controlcard,"NDIST",ndist_,0)
2724       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2725       if (restr_type.eq.10)
2726      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2727       if (restr_type.eq.12)
2728      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2729       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2730       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2731       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2732       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2733       normalize = index(controlcard,"NORMALIZE").gt.0
2734       write (iout,*) "WBOLTZD",wboltzd
2735       write (iout,*) "SCAL_PEAK",scal_peak
2736       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2737       write (iout,*) "IFRAG"
2738       do i=1,nfrag_
2739         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2740       enddo
2741       write (iout,*) "IPAIR"
2742       do i=1,npair_
2743         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2744       enddo
2745       if (nfrag_.gt.0) write (iout,*)
2746      &   "Distance restraints as generated from reference structure"
2747       do i=1,nfrag_
2748         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2749         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2750      &    ifrag_(2,i)=nstart_sup+nsup-1
2751 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2752 c        call flush(iout)
2753         if (wfrag_(i).eq.0.0d0) cycle
2754         do j=ifrag_(1,i),ifrag_(2,i)-1
2755           do k=j+1,ifrag_(2,i)
2756 c            write (iout,*) "j",j," k",k
2757             ddjk=dist(j,k)
2758             if (restr_type.eq.1) then
2759               nhpb=nhpb+1
2760               irestr_type(nhpb)=1
2761               ihpb(nhpb)=j
2762               jhpb(nhpb)=k
2763               dhpb(nhpb)=ddjk
2764               forcon(nhpb)=wfrag_(i)
2765             else if (constr_dist.eq.2) then
2766               if (ddjk.le.dist_cut) then
2767                 nhpb=nhpb+1
2768                 irestr_type(nhpb)=1
2769                 ihpb(nhpb)=j
2770                 jhpb(nhpb)=k
2771                 dhpb(nhpb)=ddjk
2772                 forcon(nhpb)=wfrag_(i)
2773               endif
2774             else if (restr_type.eq.3) then
2775               nhpb=nhpb+1
2776               irestr_type(nhpb)=1
2777               ihpb(nhpb)=j
2778               jhpb(nhpb)=k
2779               dhpb(nhpb)=ddjk
2780               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2781             endif
2782 #ifdef MPI
2783             if (.not.out1file .or. me.eq.king)
2784      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2785      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2786 #else
2787             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2788      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2789 #endif
2790           enddo
2791         enddo
2792       enddo
2793       do i=1,npair_
2794         if (wpair_(i).eq.0.0d0) cycle
2795         ii = ipair_(1,i)
2796         jj = ipair_(2,i)
2797         if (ii.gt.jj) then
2798           itemp=ii
2799           ii=jj
2800           jj=itemp
2801         endif
2802         do j=ifrag_(1,ii),ifrag_(2,ii)
2803           do k=ifrag_(1,jj),ifrag_(2,jj)
2804             ddjk=dist(j,k)
2805             if (restr_type.eq.1) then
2806               nhpb=nhpb+1
2807               irestr_type(nhpb)=1
2808               ihpb(nhpb)=j
2809               jhpb(nhpb)=k
2810               dhpb(nhpb)=ddjk
2811               forcon(nhpb)=wpair_(i)
2812             else if (constr_dist.eq.2) then
2813               if (ddjk.le.dist_cut) then
2814                 nhpb=nhpb+1
2815                 irestr_type(nhpb)=1
2816                 ihpb(nhpb)=j
2817                 jhpb(nhpb)=k
2818                 dhpb(nhpb)=ddjk
2819                 forcon(nhpb)=wpair_(i)
2820               endif
2821             else if (restr_type.eq.3) then
2822               nhpb=nhpb+1
2823               irestr_type(nhpb)=1
2824               ihpb(nhpb)=j
2825               jhpb(nhpb)=k
2826               dhpb(nhpb)=ddjk
2827               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2828             endif
2829 #ifdef MPI
2830             if (.not.out1file .or. me.eq.king)
2831      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2832      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2833 #else
2834             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2835      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2836 #endif
2837           enddo
2838         enddo
2839       enddo
2840
2841 c      print *,ndist_
2842       write (iout,*) "Distance restraints as read from input"
2843       do i=1,ndist_
2844         if (restr_type.eq.12) then
2845           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2846      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2847      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2848      &    fordepth_peak(nhpb_peak+1),npeak
2849 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2850 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2851 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2852 c     &    fordepth_peak(nhpb_peak+1),npeak
2853           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2854      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2855           nhpb_peak=nhpb_peak+1
2856           irestr_type_peak(nhpb_peak)=12
2857           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2858           ipeak(2,npeak)=i
2859 #ifdef MPI
2860           if (.not.out1file .or. me.eq.king)
2861      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2862      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2863      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2864      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2865      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2866 #else
2867           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2868      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2869      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2870      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2871      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2872 #endif
2873           if (ibecarb_peak(nhpb_peak).gt.0) then
2874             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2875             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2876           endif
2877         else if (restr_type.eq.11) then
2878           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2879      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2880 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2881           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2882           nhpb=nhpb+1
2883           irestr_type(nhpb)=11
2884 #ifdef MPI
2885           if (.not.out1file .or. me.eq.king)
2886      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2887      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2888      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2889 #else
2890           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2891      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2892      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2893 #endif
2894           if (ibecarb(nhpb).gt.0) then
2895             ihpb(nhpb)=ihpb(nhpb)+nres
2896             jhpb(nhpb)=jhpb(nhpb)+nres
2897           endif
2898         else if (restr_type.eq.10) then
2899 c Cross-lonk Markov-like potential
2900           call card_concat(controlcard)
2901           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2902           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2903           ibecarb(nhpb+1)=0
2904           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2905           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2906           if (index(controlcard,"ZL").gt.0) then
2907             link_type=1
2908           else if (index(controlcard,"ADH").gt.0) then
2909             link_type=2
2910           else if (index(controlcard,"PDH").gt.0) then
2911             link_type=3
2912           else if (index(controlcard,"DSS").gt.0) then
2913             link_type=4
2914           else
2915             link_type=0
2916           endif
2917           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2918      &       xlink(1,link_type))
2919           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2920      &       xlink(2,link_type))
2921           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2922      &       xlink(3,link_type))
2923           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2924      &       xlink(4,link_type))
2925           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2926 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2927 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2928           if (forcon(nhpb+1).le.0.0d0 .or.
2929      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2930           nhpb=nhpb+1
2931           irestr_type(nhpb)=10
2932           if (ibecarb(nhpb).gt.0) then
2933             ihpb(nhpb)=ihpb(nhpb)+nres
2934             jhpb(nhpb)=jhpb(nhpb)+nres
2935           endif
2936 #ifdef MPI
2937           if (.not.out1file .or. me.eq.king)
2938      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2939      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2940      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2941      &     irestr_type(nhpb)
2942 #else
2943           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2944      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2945      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2946      &     irestr_type(nhpb)
2947 #endif
2948         else
2949 C        print *,"in else"
2950           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2951      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2952           if (forcon(nhpb+1).gt.0.0d0) then
2953           nhpb=nhpb+1
2954           if (dhpb1(nhpb).eq.0.0d0) then
2955             irestr_type(nhpb)=1
2956           else
2957             irestr_type(nhpb)=2
2958           endif
2959           if (ibecarb(nhpb).gt.0) then
2960             ihpb(nhpb)=ihpb(nhpb)+nres
2961             jhpb(nhpb)=jhpb(nhpb)+nres
2962           endif
2963           if (dhpb(nhpb).eq.0.0d0)
2964      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2965         endif
2966 #ifdef MPI
2967           if (.not.out1file .or. me.eq.king)
2968      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2969      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2970 #else
2971           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2972      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2973 #endif
2974         endif
2975 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2976 C        if (forcon(nhpb+1).gt.0.0d0) then
2977 C          nhpb=nhpb+1
2978 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2979       enddo
2980
2981       ENDDO ! next
2982
2983       fordepthmax=0.0d0
2984       if (normalize) then
2985         do i=nss+1,nhpb
2986           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2987      &      fordepthmax=fordepth(i)
2988         enddo
2989         do i=nss+1,nhpb
2990           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2991         enddo
2992       endif
2993       return
2994       end
2995 c-------------------------------------------------------------------------------
2996 #ifdef WINIFL
2997       subroutine flush(iu)
2998       return
2999       end
3000 #endif
3001 #ifdef AIX
3002       subroutine flush(iu)
3003       call flush_(iu)
3004       return
3005       end
3006 #endif
3007 c------------------------------------------------------------------------------
3008       subroutine copy_to_tmp(source)
3009       include "DIMENSIONS"
3010       include "COMMON.IOUNITS"
3011       character*(*) source
3012       character* 256 tmpfile
3013       integer ilen
3014       external ilen
3015       logical ex
3016       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3017       inquire(file=tmpfile,exist=ex)
3018       if (ex) then
3019         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3020      &   " to temporary directory..."
3021         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3022         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3023       endif
3024       return
3025       end
3026 c------------------------------------------------------------------------------
3027       subroutine move_from_tmp(source)
3028       include "DIMENSIONS"
3029       include "COMMON.IOUNITS"
3030       character*(*) source
3031       integer ilen
3032       external ilen
3033       write (*,*) "Moving ",source(:ilen(source)),
3034      & " from temporary directory to working directory"
3035       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3036       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3037       return
3038       end
3039 c------------------------------------------------------------------------------
3040       subroutine random_init(seed)
3041 C
3042 C Initialize random number generator
3043 C
3044       implicit real*8 (a-h,o-z)
3045       include 'DIMENSIONS'
3046 #ifdef MPI
3047       include 'mpif.h'
3048       logical OKRandom, prng_restart
3049       real*8  r1
3050       integer iseed_array(4)
3051 #endif
3052       include 'COMMON.IOUNITS'
3053       include 'COMMON.TIME1'
3054       include 'COMMON.THREAD'
3055       include 'COMMON.SBRIDGE'
3056       include 'COMMON.CONTROL'
3057       include 'COMMON.MCM'
3058       include 'COMMON.MAP'
3059       include 'COMMON.HEADER'
3060       include 'COMMON.CSA'
3061       include 'COMMON.CHAIN'
3062       include 'COMMON.MUCA'
3063       include 'COMMON.MD'
3064       include 'COMMON.FFIELD'
3065       include 'COMMON.SETUP'
3066       iseed=-dint(dabs(seed))
3067       if (iseed.eq.0) then
3068         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3069      &    'Random seed undefined. The program will stop.'
3070         write (*,'(/80(1h*)/20x,a/80(1h*))')
3071      &    'Random seed undefined. The program will stop.'
3072 #ifdef MPI
3073         call mpi_finalize(mpi_comm_world,ierr)
3074 #endif
3075         stop 'Bad random seed.'
3076       endif
3077 #ifdef MPI
3078       if (fg_rank.eq.0) then
3079       seed=seed*(me+1)+1
3080 #ifdef AMD64
3081       if(me.eq.king)
3082      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3083       OKRandom = prng_restart(me,iseed)
3084 #else
3085       do i=1,4
3086        tmp=65536.0d0**(4-i)
3087        iseed_array(i) = dint(seed/tmp)
3088        seed=seed-iseed_array(i)*tmp
3089       enddo
3090       if(me.eq.king)
3091      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3092      &                 (iseed_array(i),i=1,4)
3093       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3094      &                 (iseed_array(i),i=1,4)
3095       OKRandom = prng_restart(me,iseed_array)
3096 #endif
3097       if (OKRandom) then
3098 c        r1 = prng_next(me)
3099         r1=ran_number(0.0D0,1.0D0)
3100         if(me.eq.king)
3101      &   write (iout,*) 'ran_num',r1
3102         if (r1.lt.0.0d0) OKRandom=.false.
3103       endif
3104       if (.not.OKRandom) then
3105         write (iout,*) 'PRNG IS NOT WORKING!!!'
3106         print *,'PRNG IS NOT WORKING!!!'
3107         if (me.eq.0) then
3108          call flush(iout)
3109          call mpi_abort(mpi_comm_world,error_msg,ierr)
3110          stop
3111         else
3112          write (iout,*) 'too many processors for parallel prng'
3113          write (*,*) 'too many processors for parallel prng'
3114          call flush(iout)
3115          stop
3116         endif
3117       endif
3118       endif
3119 #else
3120       call vrndst(iseed)
3121       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3122 #endif
3123       return
3124       end