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