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