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