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

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