source/unres/src-HCD-5D/readrtns_CSA.F

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