2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read MCM control parameters if required
19 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 csa if (modecalc.eq.8) then
32 csa inquire (file="fort.40",exist=file_exist)
33 csa if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
48 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
76 csa include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
180 write(*,*) "CSA not supported in this version"
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
219 if(me.eq.king.or..not.out1file)
220 & write (iout,'(2a)') diagmeth(kdiag),
221 & ' routine used to diagonalize matrices.'
224 c--------------------------------------------------------------------------
225 subroutine read_REMDpar
229 implicit real*8 (a-h,o-z)
231 include 'COMMON.IOUNITS'
232 include 'COMMON.TIME1'
235 include 'COMMON.LANGEVIN'
237 include 'COMMON.LANGEVIN.lang0'
239 include 'COMMON.INTERACT'
240 include 'COMMON.NAMES'
242 include 'COMMON.REMD'
243 include 'COMMON.CONTROL'
244 include 'COMMON.SETUP'
246 character*320 controlcard
247 character*3200 controlcard1
248 integer iremd_m_total
250 if(me.eq.king.or..not.out1file)
251 & write (iout,*) "REMD setup"
253 call card_concat(controlcard)
254 call readi(controlcard,"NREP",nrep,3)
255 call readi(controlcard,"NSTEX",nstex,1000)
256 call reada(controlcard,"RETMIN",retmin,10.0d0)
257 call reada(controlcard,"RETMAX",retmax,1000.0d0)
258 mremdsync=(index(controlcard,'SYNC').gt.0)
259 call readi(controlcard,"NSYN",i_sync_step,100)
260 restart1file=(index(controlcard,'REST1FILE').gt.0)
261 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
262 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
263 if(max_cache_traj_use.gt.max_cache_traj)
264 & max_cache_traj_use=max_cache_traj
265 if(me.eq.king.or..not.out1file) then
266 cd if (traj1file) then
267 crc caching is in testing - NTWX is not ignored
268 cd write (iout,*) "NTWX value is ignored"
269 cd write (iout,*) " trajectory is stored to one file by master"
270 cd write (iout,*) " before exchange at NSTEX intervals"
272 write (iout,*) "NREP= ",nrep
273 write (iout,*) "NSTEX= ",nstex
274 write (iout,*) "SYNC= ",mremdsync
275 write (iout,*) "NSYN= ",i_sync_step
276 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
279 t_exchange_only=(index(controlcard,'TONLY').gt.0)
280 call readi(controlcard,"HREMD",hremd,0)
281 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
282 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
284 if(usampl.and.hremd.gt.0) then
286 & "========== ERROR: USAMPL and HREMD cannot be used together"
288 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
295 if (index(controlcard,'TLIST').gt.0) then
297 call card_concat(controlcard1)
298 read(controlcard1,*) (remd_t(i),i=1,nrep)
299 if(me.eq.king.or..not.out1file)
300 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
303 if (index(controlcard,'MLIST').gt.0) then
305 call card_concat(controlcard1)
306 read(controlcard1,*) (remd_m(i),i=1,nrep)
307 if(me.eq.king.or..not.out1file) then
308 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
311 iremd_m_total=iremd_m_total+remd_m(i)
314 write (iout,*) 'Total number of replicas ',
315 & iremd_m_total*hremd
317 write (iout,*) 'Total number of replicas ',iremd_m_total
321 if(me.eq.king.or..not.out1file)
322 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
325 c--------------------------------------------------------------------------
326 subroutine read_MDpar
330 implicit real*8 (a-h,o-z)
332 include 'COMMON.IOUNITS'
333 include 'COMMON.TIME1'
336 include 'COMMON.LANGEVIN'
338 include 'COMMON.LANGEVIN.lang0'
340 include 'COMMON.INTERACT'
341 include 'COMMON.NAMES'
343 include 'COMMON.SETUP'
344 include 'COMMON.CONTROL'
345 include 'COMMON.SPLITELE'
347 character*320 controlcard
349 call card_concat(controlcard)
350 call readi(controlcard,"NSTEP",n_timestep,1000000)
351 call readi(controlcard,"NTWE",ntwe,100)
352 call readi(controlcard,"NTWX",ntwx,1000)
353 call reada(controlcard,"DT",d_time,1.0d-1)
354 call reada(controlcard,"DVMAX",dvmax,2.0d1)
355 call reada(controlcard,"DAMAX",damax,1.0d1)
356 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
357 call readi(controlcard,"LANG",lang,0)
358 RESPA = index(controlcard,"RESPA") .gt. 0
359 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
360 ntime_split0=ntime_split
361 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
362 ntime_split0=ntime_split
363 call reada(controlcard,"R_CUT",r_cut,2.0d0)
364 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
365 rest = index(controlcard,"REST").gt.0
366 tbf = index(controlcard,"TBF").gt.0
367 call readi(controlcard,"HMC",hmc,0)
368 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
369 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
370 tnh = index(controlcard,"NOSEHOOVER96").gt.0
371 if (RESPA.and.tnh)then
372 xiresp = index(controlcard,"XIRESP").gt.0
374 call reada(controlcard,"Q_NP",Q_np,0.1d0)
375 usampl = index(controlcard,"USAMPL").gt.0
377 mdpdb = index(controlcard,"MDPDB").gt.0
378 call reada(controlcard,"T_BATH",t_bath,300.0d0)
379 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
380 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
381 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
382 if (count_reset_moment.eq.0) count_reset_moment=1000000000
383 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
384 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
385 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
386 if (count_reset_vel.eq.0) count_reset_vel=1000000000
387 large = index(controlcard,"LARGE").gt.0
388 print_compon = index(controlcard,"PRINT_COMPON").gt.0
389 rattle = index(controlcard,"RATTLE").gt.0
390 c if performing umbrella sampling, fragments constrained are read from the fragment file
396 if(me.eq.king.or..not.out1file) then
398 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
400 write (iout,'(a)') "The units are:"
401 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
402 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
403 & " acceleration: angstrom/(48.9 fs)**2"
404 write (iout,'(a)') "energy: kcal/mol, temperature: K"
406 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
407 write (iout,'(a60,f10.5,a)')
408 & "Initial time step of numerical integration:",d_time,
410 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
412 write (iout,'(2a,i4,a)')
413 & "A-MTS algorithm used; initial time step for fast-varying",
414 & " short-range forces split into",ntime_split," steps."
415 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
416 & r_cut," lambda",rlamb
418 write (iout,'(2a,f10.5)')
419 & "Maximum acceleration threshold to reduce the time step",
420 & "/increase split number:",damax
421 write (iout,'(2a,f10.5)')
422 & "Maximum predicted energy drift to reduce the timestep",
423 & "/increase split number:",edriftmax
424 write (iout,'(a60,f10.5)')
425 & "Maximum velocity threshold to reduce velocities:",dvmax
426 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
427 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
428 if (rattle) write (iout,'(a60)')
429 & "Rattle algorithm used to constrain the virtual bonds"
433 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
434 call reada(controlcard,"RWAT",rwat,1.4d0)
435 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
436 surfarea=index(controlcard,"SURFAREA").gt.0
437 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
438 if(me.eq.king.or..not.out1file)then
439 write (iout,'(/a,$)') "Langevin dynamics calculation"
442 & " with direct integration of Langevin equations"
443 else if (lang.eq.2) then
444 write (iout,'(a/)') " with TINKER stochasic MD integrator"
445 else if (lang.eq.3) then
446 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
447 else if (lang.eq.4) then
448 write (iout,'(a/)') " in overdamped mode"
450 write (iout,'(//a,i5)')
451 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
454 write (iout,'(a60,f10.5)') "Temperature:",t_bath
455 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
456 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
457 write (iout,'(a60,f10.5)')
458 & "Scaling factor of the friction forces:",scal_fric
459 if (surfarea) write (iout,'(2a,i10,a)')
460 & "Friction coefficients will be scaled by solvent-accessible",
461 & " surface area every",reset_fricmat," steps."
463 c Calculate friction coefficients and bounds of stochastic forces
464 eta=6*pi*cPoise*etawat
465 if(me.eq.king.or..not.out1file)
466 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
468 gamp=scal_fric*(pstok+rwat)*eta
469 stdfp=dsqrt(2*Rb*t_bath/d_time)
471 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
472 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
474 if(me.eq.king.or..not.out1file)then
475 write (iout,'(/2a/)')
476 & "Radii of site types and friction coefficients and std's of",
477 & " stochastic forces of fully exposed sites"
478 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
480 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
481 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
485 if(me.eq.king.or..not.out1file)then
486 write (iout,'(a)') "Berendsen bath calculation"
487 write (iout,'(a60,f10.5)') "Temperature:",t_bath
488 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
490 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
491 & count_reset_moment," steps"
493 & write (iout,'(a,i10,a)')
494 & "Velocities will be reset at random every",count_reset_vel,
497 else if (tnp .or. tnp1 .or. tnh) then
498 if (tnp .or. tnp1) then
499 write (iout,'(a)') "Nose-Poincare bath calculation"
500 if (tnp) write (iout,'(a)')
501 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
502 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
504 write (iout,'(a)') "Nose-Hoover bath calculation"
505 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
515 WDTI(i) = 1.0*d_time/nresn
522 write (iout,'(a)') "NVT-XI-RESPA algorithm"
524 write (iout,'(a)') "NVT-XO-RESPA algorithm"
527 WDTIi(i) = 1.0*d_time/nresn/ntime_split
535 write (iout,'(a60,f10.5)') "Temperature:",t_bath
536 write (iout,'(a60,f10.5)') "Q =",Q_np
538 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
539 & count_reset_moment," steps"
541 & write (iout,'(a,i10,a)')
542 & "Velocities will be reset at random every",count_reset_vel,
545 else if (hmc.gt.0) then
546 write (iout,'(a)') "Hybrid Monte Carlo calculation"
547 write (iout,'(a60,f10.5)') "Temperature:",t_bath
548 write (iout,'(a60,i10)')
549 & "Number of MD steps between Metropolis tests:",hmc
552 if(me.eq.king.or..not.out1file)
553 & write (iout,'(a31)') "Microcanonical mode calculation"
555 if(me.eq.king.or..not.out1file)then
556 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
558 write(iout,*) "MD running with constraints."
559 write(iout,*) "Equilibration time ", eq_time, " mtus."
560 write(iout,*) "Constraining ", nfrag," fragments."
561 write(iout,*) "Length of each fragment, weight and q0:"
563 write (iout,*) "Set of restraints #",iset
565 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
566 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
568 write(iout,*) "constraints between ", npair, "fragments."
569 write(iout,*) "constraint pairs, weights and q0:"
571 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
572 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
574 write(iout,*) "angle constraints within ", nfrag_back,
575 & "backbone fragments."
576 write(iout,*) "fragment, weights:"
578 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
579 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
580 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
583 iset=mod(kolor,nset)+1
586 if(me.eq.king.or..not.out1file)
587 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
590 c------------------------------------------------------------------------------
593 C Read molecular data.
595 implicit real*8 (a-h,o-z)
601 include 'COMMON.IOUNITS'
604 include 'COMMON.INTERACT'
605 include 'COMMON.LOCAL'
606 include 'COMMON.NAMES'
607 include 'COMMON.CHAIN'
608 include 'COMMON.FFIELD'
609 include 'COMMON.SBRIDGE'
610 include 'COMMON.HEADER'
611 include 'COMMON.CONTROL'
612 include 'COMMON.DBASE'
613 include 'COMMON.THREAD'
614 include 'COMMON.CONTACTS'
615 include 'COMMON.TORCNSTR'
616 include 'COMMON.TIME1'
617 include 'COMMON.BOUNDS'
619 include 'COMMON.REMD'
620 include 'COMMON.SETUP'
621 character*4 sequence(maxres)
623 double precision x(maxvar)
624 character*256 pdbfile
625 character*320 weightcard
626 character*80 weightcard_t,ucase
627 dimension itype_pdb(maxres)
628 common /pizda/ itype_pdb
629 logical seq_comp,fail
630 double precision energia(0:n_ene)
636 C Read weights of the subsequent energy terms.
649 if(me.eq.king.or..not.out1file) then
650 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
651 write (iout,*) 'Current weights for processor ',
652 & me,' set ',i2set(me)
656 call card_concat(weightcard)
657 call reada(weightcard,'WLONG',wlong,1.0D0)
658 call reada(weightcard,'WSC',wsc,wlong)
659 call reada(weightcard,'WSCP',wscp,wlong)
660 call reada(weightcard,'WELEC',welec,1.0D0)
661 call reada(weightcard,'WVDWPP',wvdwpp,welec)
662 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
663 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
664 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
665 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
666 call reada(weightcard,'WTURN3',wturn3,1.0D0)
667 call reada(weightcard,'WTURN4',wturn4,1.0D0)
668 call reada(weightcard,'WTURN6',wturn6,1.0D0)
669 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
670 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
671 call reada(weightcard,'WBOND',wbond,1.0D0)
672 call reada(weightcard,'WTOR',wtor,1.0D0)
673 call reada(weightcard,'WTORD',wtor_d,1.0D0)
674 call reada(weightcard,'WANG',wang,1.0D0)
675 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
676 call reada(weightcard,'SCAL14',scal14,0.4D0)
677 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
678 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
679 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
680 call reada(weightcard,'TEMP0',temp0,300.0d0)
681 if (index(weightcard,'SOFT').gt.0) ipot=6
682 C 12/1/95 Added weight for the multi-body term WCORR
683 call reada(weightcard,'WCORRH',wcorr,1.0D0)
684 if (wcorr4.gt.0.0d0) wcorr=wcorr4
692 hweights(i,7)=wel_loc
695 hweights(i,10)=wturn6
697 hweights(i,12)=wscloc
699 hweights(i,14)=wtor_d
700 hweights(i,15)=wstrain
701 hweights(i,16)=wvdwpp
703 hweights(i,18)=scal14
704 hweights(i,21)=wsccor
709 weights(i)=hweights(i2set(me),i)
733 call card_concat(weightcard)
734 call reada(weightcard,'WLONG',wlong,1.0D0)
735 call reada(weightcard,'WSC',wsc,wlong)
736 call reada(weightcard,'WSCP',wscp,wlong)
737 call reada(weightcard,'WELEC',welec,1.0D0)
738 call reada(weightcard,'WVDWPP',wvdwpp,welec)
739 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
740 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
741 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
742 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
743 call reada(weightcard,'WTURN3',wturn3,1.0D0)
744 call reada(weightcard,'WTURN4',wturn4,1.0D0)
745 call reada(weightcard,'WTURN6',wturn6,1.0D0)
746 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
747 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
748 call reada(weightcard,'WBOND',wbond,1.0D0)
749 call reada(weightcard,'WTOR',wtor,1.0D0)
750 call reada(weightcard,'WTORD',wtor_d,1.0D0)
751 call reada(weightcard,'WANG',wang,1.0D0)
752 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
753 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
754 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
755 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
756 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
757 call reada(weightcard,'SCAL14',scal14,0.4D0)
758 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
759 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
760 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
761 call reada(weightcard,'TEMP0',temp0,300.0d0)
762 if (index(weightcard,'SOFT').gt.0) ipot=6
763 C 12/1/95 Added weight for the multi-body term WCORR
764 call reada(weightcard,'WCORRH',wcorr,1.0D0)
765 if (wcorr4.gt.0.0d0) wcorr=wcorr4
786 weights(25)=wdfa_dist
789 weights(28)=wdfa_beta
791 if(me.eq.king.or..not.out1file)
792 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
793 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
795 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
796 10 format (/'Energy-term weights (unscaled):'//
797 & 'WSCC= ',f10.6,' (SC-SC)'/
798 & 'WSCP= ',f10.6,' (SC-p)'/
799 & 'WELEC= ',f10.6,' (p-p electr)'/
800 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
801 & 'WBOND= ',f10.6,' (stretching)'/
802 & 'WANG= ',f10.6,' (bending)'/
803 & 'WSCLOC= ',f10.6,' (SC local)'/
804 & 'WTOR= ',f10.6,' (torsional)'/
805 & 'WTORD= ',f10.6,' (double torsional)'/
806 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
807 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
808 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
809 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
810 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
811 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
812 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
813 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
814 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
815 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
816 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
817 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
818 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
819 if(me.eq.king.or..not.out1file)then
820 if (wcorr4.gt.0.0d0) then
821 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
822 & 'between contact pairs of peptide groups'
823 write (iout,'(2(a,f5.3/))')
824 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
825 & 'Range of quenching the correlation terms:',2*delt_corr
826 else if (wcorr.gt.0.0d0) then
827 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
828 & 'between contact pairs of peptide groups'
830 write (iout,'(a,f8.3)')
831 & 'Scaling factor of 1,4 SC-p interactions:',scal14
832 write (iout,'(a,f8.3)')
833 & 'General scaling factor of SC-p interactions:',scalscp
835 r0_corr=cutoff_corr-delt_corr
837 aad(i,1)=scalscp*aad(i,1)
838 aad(i,2)=scalscp*aad(i,2)
839 bad(i,1)=scalscp*bad(i,1)
840 bad(i,2)=scalscp*bad(i,2)
842 call rescale_weights(t_bath)
843 if(me.eq.king.or..not.out1file)
844 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
845 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
847 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
848 22 format (/'Energy-term weights (scaled):'//
849 & 'WSCC= ',f10.6,' (SC-SC)'/
850 & 'WSCP= ',f10.6,' (SC-p)'/
851 & 'WELEC= ',f10.6,' (p-p electr)'/
852 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
853 & 'WBOND= ',f10.6,' (stretching)'/
854 & 'WANG= ',f10.6,' (bending)'/
855 & 'WSCLOC= ',f10.6,' (SC local)'/
856 & 'WTOR= ',f10.6,' (torsional)'/
857 & 'WTORD= ',f10.6,' (double torsional)'/
858 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
859 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
860 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
861 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
862 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
863 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
864 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
865 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
866 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
867 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
868 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
869 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
870 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
871 if(me.eq.king.or..not.out1file)
872 & write (iout,*) "Reference temperature for weights calculation:",
874 call reada(weightcard,"D0CM",d0cm,3.78d0)
875 call reada(weightcard,"AKCM",akcm,15.1d0)
876 call reada(weightcard,"AKTH",akth,11.0d0)
877 call reada(weightcard,"AKCT",akct,12.0d0)
878 call reada(weightcard,"V1SS",v1ss,-1.08d0)
879 call reada(weightcard,"V2SS",v2ss,7.61d0)
880 call reada(weightcard,"V3SS",v3ss,13.7d0)
881 call reada(weightcard,"EBR",ebr,-5.50D0)
882 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
884 dyn_ss_mask(i)=.false.
888 dyn_ssbond_ij(i,j)=1.0d300
891 call reada(weightcard,"HT",Ht,0.0D0)
893 ss_depth=ebr/wsc-0.25*eps(1,1)
894 Ht=Ht/wsc-0.25*eps(1,1)
895 akcm=akcm*wstrain/wsc
896 akth=akth*wstrain/wsc
897 akct=akct*wstrain/wsc
898 v1ss=v1ss*wstrain/wsc
899 v2ss=v2ss*wstrain/wsc
900 v3ss=v3ss*wstrain/wsc
902 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
905 if(me.eq.king.or..not.out1file) then
906 write (iout,*) "Parameters of the SS-bond potential:"
907 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
909 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
910 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
911 write (iout,*)" HT",Ht
912 print *,'indpdb=',indpdb,' pdbref=',pdbref
914 if (indpdb.gt.0 .or. pdbref) then
915 read(inp,'(a)') pdbfile
916 if(me.eq.king.or..not.out1file)
917 & write (iout,'(2a)') 'PDB data will be read from file ',
918 & pdbfile(:ilen(pdbfile))
919 open(ipdbin,file=pdbfile,status='old',err=33)
921 33 write (iout,'(a)') 'Error opening PDB file.'
924 c print *,'Begin reading pdb data'
926 c print *,'Finished reading pdb data'
927 if(me.eq.king.or..not.out1file)
928 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
929 & ' nstart_sup=',nstart_sup
931 itype_pdb(i)=itype(i)
935 nct=nstart_sup+nsup-1
936 call contact(.false.,ncont_ref,icont_ref,co)
939 C Following 2 lines for diagnostics; comment out if not needed
940 write (iout,*) "Before sideadd"
942 if(me.eq.king.or..not.out1file)
943 & write(iout,*)'Adding sidechains'
950 do while (fail.and.nsi.le.maxsi)
951 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
954 if(fail) write(iout,*)'Adding sidechain failed for res ',
955 & i,' after ',nsi,' trials'
958 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
961 C Following 2 lines for diagnostics; comment out if not needed
962 c write (iout,*) "After sideadd"
965 if (indpdb.eq.0) then
966 C Read sequence if not taken from the pdb file.
968 c print *,'nres=',nres
969 if (iscode.gt.0) then
970 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
972 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
974 C Convert sequence to numeric code
976 itype(i)=rescode(i,sequence(i),iscode)
978 C Assign initial virtual bond lengths
984 vbld(i+nres)=dsc(itype(i))
985 vbld_inv(i+nres)=dsc_inv(itype(i))
986 c write (iout,*) "i",i," itype",itype(i),
987 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
991 c print '(20i4)',(itype(i),i=1,nres)
994 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
996 if (itype(i).eq.21) then
1000 else if (itype(i+1).ne.20) then
1002 else if (itype(i).ne.20) then
1009 if(me.eq.king.or..not.out1file)then
1010 write (iout,*) "ITEL"
1012 write (iout,*) i,itype(i),itel(i)
1014 print *,'Call Read_Bridge.'
1017 C 8/13/98 Set limits to generating the dihedral angles
1022 read (inp,*) ndih_constr
1023 if (ndih_constr.gt.0) then
1025 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1026 if(me.eq.king.or..not.out1file)then
1028 & 'There are',ndih_constr,' constraints on phi angles.'
1030 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1034 phi0(i)=deg2rad*phi0(i)
1035 drange(i)=deg2rad*drange(i)
1037 if(me.eq.king.or..not.out1file)
1038 & write (iout,*) 'FTORS',ftors
1041 phibound(1,ii) = phi0(i)-drange(i)
1042 phibound(2,ii) = phi0(i)+drange(i)
1047 if (me.eq.king) then
1049 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1051 write (iout,'(a3,i5,2f10.1)')
1052 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1058 cd print *,'NNT=',NNT,' NCT=',NCT
1059 if (itype(1).eq.21) nnt=2
1060 if (itype(nres).eq.21) nct=nct-1
1062 C Bartek:READ init_vars
1063 C Initialize variables!
1064 C Juyong:READ read_info
1065 C READ fragment information!!
1066 C both routines should be in dfa.F file!!
1068 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1069 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1071 print*, 'init_dfa_vars finished!'
1073 print*, 'read_dfa_info finished!'
1077 if(me.eq.king.or..not.out1file)
1078 & write (iout,'(a,i3)') 'nsup=',nsup
1080 if (nsup.le.(nct-nnt+1)) then
1081 do i=0,nct-nnt+1-nsup
1082 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1088 & 'Error - sequences to be superposed do not match.'
1091 do i=0,nsup-(nct-nnt+1)
1092 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1094 nstart_sup=nstart_sup+i
1100 & 'Error - sequences to be superposed do not match.'
1103 if (nsup.eq.0) nsup=nct-nnt
1104 if (nstart_sup.eq.0) nstart_sup=nnt
1105 if (nstart_seq.eq.0) nstart_seq=nnt
1106 if(me.eq.king.or..not.out1file)
1107 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1108 & ' nstart_seq=',nstart_seq
1110 c--- Zscore rms -------
1111 if (nz_start.eq.0) nz_start=nnt
1112 if (nz_end.eq.0 .and. nsup.gt.0) then
1114 else if (nz_end.eq.0) then
1117 if(me.eq.king.or..not.out1file)then
1118 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1119 write (iout,*) 'IZ_SC=',iz_sc
1121 c----------------------
1124 if (.not.pdbref) then
1125 call read_angles(inp,*38)
1127 38 write (iout,'(a)') 'Error reading reference structure.'
1129 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1130 stop 'Error reading reference structure'
1134 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1143 call contact(.true.,ncont_ref,icont_ref,co)
1145 if(me.eq.king.or..not.out1file)
1146 & write (iout,*) 'Contact order:',co
1148 if(me.eq.king.or..not.out1file)
1149 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1152 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1154 if(me.eq.king.or..not.out1file)
1155 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1156 & icont_ref(1,i),' ',
1157 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1161 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1162 if (constr_dist.gt.0) then
1163 call read_dist_constr
1167 if (constr_homology.gt.0) then
1168 call read_constr_homology
1174 if (nhpb.gt.0) call hpb_partition
1175 c write (iout,*) "After read_dist_constr nhpb",nhpb
1177 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1178 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1179 & modecalc.ne.10) then
1180 C If input structure hasn't been supplied from the PDB file read or generate
1182 if (iranconf.eq.0 .and. .not. extconf) then
1183 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1184 & write (iout,'(a)') 'Initial geometry will be read in.'
1186 read(inp,'(8f10.5)',end=36,err=36)
1187 & ((c(l,k),l=1,3),k=1,nres),
1188 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1189 call int_from_cart1(.false.)
1192 dc(j,i)=c(j,i+1)-c(j,i)
1193 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1197 if (itype(i).ne.10) then
1199 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1200 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1206 call read_angles(inp,*36)
1209 36 write (iout,'(a)') 'Error reading angle file.'
1211 call mpi_finalize( MPI_COMM_WORLD,IERR )
1213 stop 'Error reading angle file.'
1215 else if (extconf) then
1216 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1217 & write (iout,'(a)') 'Extended chain initial geometry.'
1219 theta(i)=90d0*deg2rad
1222 phi(i)=180d0*deg2rad
1225 alph(i)=110d0*deg2rad
1228 omeg(i)=-120d0*deg2rad
1231 if(me.eq.king.or..not.out1file)
1232 & write (iout,'(a)') 'Random-generated initial geometry.'
1236 if (me.eq.king .or. fg_rank.eq.0 .and. (
1237 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1241 call gen_rand_conf(itmp,*30)
1243 30 write (iout,*) 'Failed to generate random conformation',
1244 & ', itrial=',itrial
1245 write (*,*) 'Processor:',me,
1246 & ' Failed to generate random conformation',
1256 write (iout,'(a,i3,a)') 'Processor:',me,
1257 & ' error in generating random conformation.'
1258 write (*,'(a,i3,a)') 'Processor:',me,
1259 & ' error in generating random conformation.'
1262 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1269 elseif (modecalc.eq.4) then
1270 read (inp,'(a)') intinname
1271 open (intin,file=intinname,status='old',err=333)
1272 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1273 & write (iout,'(a)') 'intinname',intinname
1274 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1276 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1278 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1280 stop 'Error opening angle file.'
1284 C Generate distance constraints, if the PDB structure is to be regularized.
1285 if (nthread.gt.0) then
1286 call read_threadbase
1289 if (me.eq.king .or. .not. out1file)
1291 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1292 write (iout,'(/a,i3,a)')
1293 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1294 write (iout,'(20i4)') (iss(i),i=1,ns)
1296 write(iout,*)"Running with dynamic disulfide-bond formation"
1298 write (iout,'(/a/)') 'Pre-formed links are:'
1304 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1305 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1311 if (ns.gt.0.and.dyn_ss) then
1315 forcon(i-nss)=forcon(i)
1322 dyn_ss_mask(iss(i))=.true.
1325 if (i2ndstr.gt.0) call secstrp2dihc
1326 c call geom_to_var(nvar,x)
1327 c call etotal(energia(0))
1328 c call enerprint(energia(0))
1329 c call briefout(0,etot)
1331 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1332 cd write (iout,'(a)') 'Variable list:'
1333 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1335 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1336 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1337 & 'Processor',myrank,': end reading molecular data.'
1341 c--------------------------------------------------------------------------
1342 logical function seq_comp(itypea,itypeb,length)
1344 integer length,itypea(length),itypeb(length)
1347 if (itypea(i).ne.itypeb(i)) then
1355 c-----------------------------------------------------------------------------
1356 subroutine read_bridge
1357 C Read information about disulfide bridges.
1358 implicit real*8 (a-h,o-z)
1359 include 'DIMENSIONS'
1363 include 'COMMON.IOUNITS'
1364 include 'COMMON.GEO'
1365 include 'COMMON.VAR'
1366 include 'COMMON.INTERACT'
1367 include 'COMMON.LOCAL'
1368 include 'COMMON.NAMES'
1369 include 'COMMON.CHAIN'
1370 include 'COMMON.FFIELD'
1371 include 'COMMON.SBRIDGE'
1372 include 'COMMON.HEADER'
1373 include 'COMMON.CONTROL'
1374 include 'COMMON.DBASE'
1375 include 'COMMON.THREAD'
1376 include 'COMMON.TIME1'
1377 include 'COMMON.SETUP'
1378 C Read bridging residues.
1379 read (inp,*) ns,(iss(i),i=1,ns)
1381 if(me.eq.king.or..not.out1file)
1382 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1383 C Check whether the specified bridging residues are cystines.
1385 if (itype(iss(i)).ne.1) then
1386 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1387 & 'Do you REALLY think that the residue ',
1388 & restyp(itype(iss(i))),i,
1389 & ' can form a disulfide bridge?!!!'
1390 write (*,'(2a,i3,a)')
1391 & 'Do you REALLY think that the residue ',
1392 & restyp(itype(iss(i))),i,
1393 & ' can form a disulfide bridge?!!!'
1395 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1400 C Read preformed bridges.
1402 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1404 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1407 C Check if the residues involved in bridges are in the specified list of
1408 C bridging residues.
1411 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1412 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1413 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1414 & ' contains residues present in other pairs.'
1415 write (*,'(a,i3,a)') 'Disulfide pair',i,
1416 & ' contains residues present in other pairs.'
1418 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1424 if (ihpb(i).eq.iss(j)) goto 10
1426 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1429 if (jhpb(i).eq.iss(j)) goto 20
1431 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1437 ihpb(i)=ihpb(i)+nres
1438 jhpb(i)=jhpb(i)+nres
1444 c----------------------------------------------------------------------------
1445 subroutine read_x(kanal,*)
1446 implicit real*8 (a-h,o-z)
1447 include 'DIMENSIONS'
1448 include 'COMMON.GEO'
1449 include 'COMMON.VAR'
1450 include 'COMMON.CHAIN'
1451 include 'COMMON.IOUNITS'
1452 include 'COMMON.CONTROL'
1453 include 'COMMON.LOCAL'
1454 include 'COMMON.INTERACT'
1455 c Read coordinates from input
1457 read(kanal,'(8f10.5)',end=10,err=10)
1458 & ((c(l,k),l=1,3),k=1,nres),
1459 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1462 c(j,2*nres)=c(j,nres)
1464 call int_from_cart1(.false.)
1467 dc(j,i)=c(j,i+1)-c(j,i)
1468 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1472 if (itype(i).ne.10) then
1474 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1475 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1483 c----------------------------------------------------------------------------
1484 subroutine read_threadbase
1485 implicit real*8 (a-h,o-z)
1486 include 'DIMENSIONS'
1487 include 'COMMON.IOUNITS'
1488 include 'COMMON.GEO'
1489 include 'COMMON.VAR'
1490 include 'COMMON.INTERACT'
1491 include 'COMMON.LOCAL'
1492 include 'COMMON.NAMES'
1493 include 'COMMON.CHAIN'
1494 include 'COMMON.FFIELD'
1495 include 'COMMON.SBRIDGE'
1496 include 'COMMON.HEADER'
1497 include 'COMMON.CONTROL'
1498 include 'COMMON.DBASE'
1499 include 'COMMON.THREAD'
1500 include 'COMMON.TIME1'
1501 C Read pattern database for threading.
1502 read (icbase,*) nseq
1504 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1505 & nres_base(2,i),nres_base(3,i)
1506 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1508 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1509 c & nres_base(2,i),nres_base(3,i)
1510 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1514 if (weidis.eq.0.0D0) weidis=0.1D0
1523 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1524 write (iout,'(a,i5)') 'nexcl: ',nexcl
1525 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1528 c------------------------------------------------------------------------------
1529 subroutine setup_var
1530 implicit real*8 (a-h,o-z)
1531 include 'DIMENSIONS'
1532 include 'COMMON.IOUNITS'
1533 include 'COMMON.GEO'
1534 include 'COMMON.VAR'
1535 include 'COMMON.INTERACT'
1536 include 'COMMON.LOCAL'
1537 include 'COMMON.NAMES'
1538 include 'COMMON.CHAIN'
1539 include 'COMMON.FFIELD'
1540 include 'COMMON.SBRIDGE'
1541 include 'COMMON.HEADER'
1542 include 'COMMON.CONTROL'
1543 include 'COMMON.DBASE'
1544 include 'COMMON.THREAD'
1545 include 'COMMON.TIME1'
1546 C Set up variable list.
1552 if (itype(i).ne.10) then
1554 ialph(i,1)=nvar+nside
1558 if (indphi.gt.0) then
1560 else if (indback.gt.0) then
1565 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1568 c----------------------------------------------------------------------------
1569 subroutine gen_dist_constr
1570 C Generate CA distance constraints.
1571 implicit real*8 (a-h,o-z)
1572 include 'DIMENSIONS'
1573 include 'COMMON.IOUNITS'
1574 include 'COMMON.GEO'
1575 include 'COMMON.VAR'
1576 include 'COMMON.INTERACT'
1577 include 'COMMON.LOCAL'
1578 include 'COMMON.NAMES'
1579 include 'COMMON.CHAIN'
1580 include 'COMMON.FFIELD'
1581 include 'COMMON.SBRIDGE'
1582 include 'COMMON.HEADER'
1583 include 'COMMON.CONTROL'
1584 include 'COMMON.DBASE'
1585 include 'COMMON.THREAD'
1586 include 'COMMON.TIME1'
1587 dimension itype_pdb(maxres)
1588 common /pizda/ itype_pdb
1590 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1591 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1592 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1594 do i=nstart_sup,nstart_sup+nsup-1
1595 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1596 cd & ' seq_pdb', restyp(itype_pdb(i))
1597 do j=i+2,nstart_sup+nsup-1
1599 ihpb(nhpb)=i+nstart_seq-nstart_sup
1600 jhpb(nhpb)=j+nstart_seq-nstart_sup
1602 dhpb(nhpb)=dist(i,j)
1605 cd write (iout,'(a)') 'Distance constraints:'
1610 cd if (ii.gt.nres) then
1615 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1616 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1617 cd & dhpb(i),forcon(i)
1621 c----------------------------------------------------------------------------
1623 implicit real*8 (a-h,o-z)
1624 include 'DIMENSIONS'
1625 include 'COMMON.MAP'
1626 include 'COMMON.IOUNITS'
1627 character*3 angid(4) /'THE','PHI','ALP','OME'/
1628 character*80 mapcard,ucase
1630 read (inp,'(a)') mapcard
1631 mapcard=ucase(mapcard)
1632 if (index(mapcard,'PHI').gt.0) then
1634 else if (index(mapcard,'THE').gt.0) then
1636 else if (index(mapcard,'ALP').gt.0) then
1638 else if (index(mapcard,'OME').gt.0) then
1641 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1642 stop 'Error - illegal variable spec in MAP card.'
1644 call readi (mapcard,'RES1',res1(imap),0)
1645 call readi (mapcard,'RES2',res2(imap),0)
1646 if (res1(imap).eq.0) then
1647 res1(imap)=res2(imap)
1648 else if (res2(imap).eq.0) then
1649 res2(imap)=res1(imap)
1651 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1653 & 'Error - illegal definition of variable group in MAP.'
1654 stop 'Error - illegal definition of variable group in MAP.'
1656 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1657 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1658 call readi(mapcard,'NSTEP',nstep(imap),0)
1659 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1661 & 'Illegal boundary and/or step size specification in MAP.'
1662 stop 'Illegal boundary and/or step size specification in MAP.'
1667 c----------------------------------------------------------------------------
1668 csa subroutine csaread
1669 csa implicit real*8 (a-h,o-z)
1670 csa include 'DIMENSIONS'
1671 csa include 'COMMON.IOUNITS'
1672 csa include 'COMMON.GEO'
1673 csa include 'COMMON.CSA'
1674 csa include 'COMMON.BANK'
1675 csa include 'COMMON.CONTROL'
1676 csa character*80 ucase
1677 csa character*620 mcmcard
1678 csa call card_concat(mcmcard)
1680 csa call readi(mcmcard,'NCONF',nconf,50)
1681 csa call readi(mcmcard,'NADD',nadd,0)
1682 csa call readi(mcmcard,'JSTART',jstart,1)
1683 csa call readi(mcmcard,'JEND',jend,1)
1684 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1685 csa call readi(mcmcard,'N0',n0,1)
1686 csa call readi(mcmcard,'N1',n1,6)
1687 csa call readi(mcmcard,'N2',n2,4)
1688 csa call readi(mcmcard,'N3',n3,0)
1689 csa call readi(mcmcard,'N4',n4,0)
1690 csa call readi(mcmcard,'N5',n5,0)
1691 csa call readi(mcmcard,'N6',n6,10)
1692 csa call readi(mcmcard,'N7',n7,0)
1693 csa call readi(mcmcard,'N8',n8,0)
1694 csa call readi(mcmcard,'N9',n9,0)
1695 csa call readi(mcmcard,'N14',n14,0)
1696 csa call readi(mcmcard,'N15',n15,0)
1697 csa call readi(mcmcard,'N16',n16,0)
1698 csa call readi(mcmcard,'N17',n17,0)
1699 csa call readi(mcmcard,'N18',n18,0)
1701 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1703 csa call readi(mcmcard,'NDIFF',ndiff,2)
1704 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1705 csa call readi(mcmcard,'IS1',is1,1)
1706 csa call readi(mcmcard,'IS2',is2,8)
1707 csa call readi(mcmcard,'NRAN0',nran0,4)
1708 csa call readi(mcmcard,'NRAN1',nran1,2)
1709 csa call readi(mcmcard,'IRR',irr,1)
1710 csa call readi(mcmcard,'NSEED',nseed,20)
1711 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1712 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1713 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1714 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1715 csa call readi(mcmcard,'ICMAX',icmax,3)
1716 csa call readi(mcmcard,'IRESTART',irestart,0)
1717 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1720 csa call reada(mcmcard,'DELE',dele,20.0d0)
1721 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1722 csa call readi(mcmcard,'IREF',iref,0)
1723 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1724 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1725 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1726 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1727 csa write (iout,*) "NCONF_IN",nconf_in
1730 c----------------------------------------------------------------------------
1731 cfmc subroutine mcmfread
1732 cfmc implicit real*8 (a-h,o-z)
1733 cfmc include 'DIMENSIONS'
1734 cfmc include 'COMMON.MCMF'
1735 cfmc include 'COMMON.IOUNITS'
1736 cfmc include 'COMMON.GEO'
1737 cfmc character*80 ucase
1738 cfmc character*620 mcmcard
1739 cfmc call card_concat(mcmcard)
1741 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1742 cfmc write(iout,*)'MAXRANT=',maxrant
1743 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1744 cfmc write(iout,*)'MAXFAM=',maxfam
1745 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1746 cfmc write(iout,*)'NNET1=',nnet1
1747 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1748 cfmc write(iout,*)'NNET2=',nnet2
1749 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1750 cfmc write(iout,*)'NNET3=',nnet3
1751 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1752 cfmc write(iout,*)'ILASTT=',ilastt
1753 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1754 cfmc write(iout,*)'MAXSTR=',maxstr
1755 cfmc maxstr_f=maxstr/maxfam
1756 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1757 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1758 cfmc write(iout,*)'NMCMF=',nmcmf
1759 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1760 cfmc write(iout,*)'IFOCUS=',ifocus
1761 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1762 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1763 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1764 cfmc write(iout,*)'INTPRT=',intprt
1765 cfmc call readi(mcmcard,'IPRT',iprt,100)
1766 cfmc write(iout,*)'IPRT=',iprt
1767 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1768 cfmc write(iout,*)'IMAXTR=',imaxtr
1769 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1770 cfmc write(iout,*)'MAXEVEN=',maxeven
1771 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1772 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1773 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1774 cfmc write(iout,*)'INIMIN=',inimin
1775 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1776 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1777 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1778 cfmc write(iout,*)'NTHREAD=',nthread
1779 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1780 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1781 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1782 cfmc write(iout,*)'MAXPERT=',maxpert
1783 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1784 cfmc write(iout,*)'IRMSD=',irmsd
1785 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1786 cfmc write(iout,*)'DENEMIN=',denemin
1787 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1788 cfmc write(iout,*)'RCUT1S=',rcut1s
1789 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1790 cfmc write(iout,*)'RCUT1E=',rcut1e
1791 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1792 cfmc write(iout,*)'RCUT2S=',rcut2s
1793 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1794 cfmc write(iout,*)'RCUT2E=',rcut2e
1795 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1796 cfmc write(iout,*)'DPERT1=',d_pert1
1797 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1798 cfmc write(iout,*)'DPERT1A=',d_pert1a
1799 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1800 cfmc write(iout,*)'DPERT2=',d_pert2
1801 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1802 cfmc write(iout,*)'DPERT2A=',d_pert2a
1803 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1804 cfmc write(iout,*)'DPERT2B=',d_pert2b
1805 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1806 cfmc write(iout,*)'DPERT2C=',d_pert2c
1807 cfmc d_pert1=deg2rad*d_pert1
1808 cfmc d_pert1a=deg2rad*d_pert1a
1809 cfmc d_pert2=deg2rad*d_pert2
1810 cfmc d_pert2a=deg2rad*d_pert2a
1811 cfmc d_pert2b=deg2rad*d_pert2b
1812 cfmc d_pert2c=deg2rad*d_pert2c
1813 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1814 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1815 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1816 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1817 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1818 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1819 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1820 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1821 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1822 cfmc write(iout,*)'RCUTINI=',rcutini
1823 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1824 cfmc write(iout,*)'GRAT=',grat
1825 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1826 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1830 c----------------------------------------------------------------------------
1832 implicit real*8 (a-h,o-z)
1833 include 'DIMENSIONS'
1834 include 'COMMON.MCM'
1835 include 'COMMON.MCE'
1836 include 'COMMON.IOUNITS'
1838 character*320 mcmcard
1839 call card_concat(mcmcard)
1840 call readi(mcmcard,'MAXACC',maxacc,100)
1841 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1842 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1843 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1844 call readi(mcmcard,'MAXREPM',maxrepm,200)
1845 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1846 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1847 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1848 call reada(mcmcard,'E_UP',e_up,5.0D0)
1849 call reada(mcmcard,'DELTE',delte,0.1D0)
1850 call readi(mcmcard,'NSWEEP',nsweep,5)
1851 call readi(mcmcard,'NSTEPH',nsteph,0)
1852 call readi(mcmcard,'NSTEPC',nstepc,0)
1853 call reada(mcmcard,'TMIN',tmin,298.0D0)
1854 call reada(mcmcard,'TMAX',tmax,298.0D0)
1855 call readi(mcmcard,'NWINDOW',nwindow,0)
1856 call readi(mcmcard,'PRINT_MC',print_mc,0)
1857 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1858 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1859 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1860 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1861 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1862 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1863 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1864 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1865 if (nwindow.gt.0) then
1866 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1868 winlen(i)=winend(i)-winstart(i)+1
1871 if (tmax.lt.tmin) tmax=tmin
1872 if (tmax.eq.tmin) then
1876 if (nstepc.gt.0 .and. nsteph.gt.0) then
1877 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1878 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1880 C Probabilities of different move types
1881 sumpro_type(0)=0.0D0
1882 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1883 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1884 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1885 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1886 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1887 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1888 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1890 print *,'i',i,' sumprotype',sumpro_type(i)
1891 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1892 print *,'i',i,' sumprotype',sumpro_type(i)
1896 c----------------------------------------------------------------------------
1897 subroutine read_minim
1898 implicit real*8 (a-h,o-z)
1899 include 'DIMENSIONS'
1900 include 'COMMON.MINIM'
1901 include 'COMMON.IOUNITS'
1903 character*320 minimcard
1904 call card_concat(minimcard)
1905 call readi(minimcard,'MAXMIN',maxmin,2000)
1906 call readi(minimcard,'MAXFUN',maxfun,5000)
1907 call readi(minimcard,'MINMIN',minmin,maxmin)
1908 call readi(minimcard,'MINFUN',minfun,maxmin)
1909 call reada(minimcard,'TOLF',tolf,1.0D-2)
1910 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1911 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1912 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1913 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1914 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1915 & 'Options in energy minimization:'
1916 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1917 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1918 & 'MinMin:',MinMin,' MinFun:',MinFun,
1919 & ' TolF:',TolF,' RTolF:',RTolF
1922 c----------------------------------------------------------------------------
1923 subroutine read_angles(kanal,*)
1924 implicit real*8 (a-h,o-z)
1925 include 'DIMENSIONS'
1926 include 'COMMON.GEO'
1927 include 'COMMON.VAR'
1928 include 'COMMON.CHAIN'
1929 include 'COMMON.IOUNITS'
1930 include 'COMMON.CONTROL'
1931 c Read angles from input
1933 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1934 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1935 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1936 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1939 c 9/7/01 avoid 180 deg valence angle
1940 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1942 theta(i)=deg2rad*theta(i)
1943 phi(i)=deg2rad*phi(i)
1944 alph(i)=deg2rad*alph(i)
1945 omeg(i)=deg2rad*omeg(i)
1950 c----------------------------------------------------------------------------
1951 subroutine reada(rekord,lancuch,wartosc,default)
1953 character*(*) rekord,lancuch
1954 double precision wartosc,default
1957 iread=index(rekord,lancuch)
1958 if (iread.eq.0) then
1962 iread=iread+ilen(lancuch)+1
1963 read (rekord(iread:),*,err=10,end=10) wartosc
1968 c----------------------------------------------------------------------------
1969 subroutine readi(rekord,lancuch,wartosc,default)
1971 character*(*) rekord,lancuch
1972 integer wartosc,default
1975 iread=index(rekord,lancuch)
1976 if (iread.eq.0) then
1980 iread=iread+ilen(lancuch)+1
1981 read (rekord(iread:),*,err=10,end=10) wartosc
1986 c----------------------------------------------------------------------------
1987 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1990 integer tablica(dim),default
1991 character*(*) rekord,lancuch
1998 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1999 if (iread.eq.0) return
2000 iread=iread+ilen(lancuch)+1
2001 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2004 c----------------------------------------------------------------------------
2005 subroutine multreada(rekord,lancuch,tablica,dim,default)
2008 double precision tablica(dim),default
2009 character*(*) rekord,lancuch
2016 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2017 if (iread.eq.0) return
2018 iread=iread+ilen(lancuch)+1
2019 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2022 c----------------------------------------------------------------------------
2023 subroutine openunits
2024 implicit real*8 (a-h,o-z)
2025 include 'DIMENSIONS'
2028 character*16 form,nodename
2031 include 'COMMON.SETUP'
2032 include 'COMMON.IOUNITS'
2034 include 'COMMON.CONTROL'
2035 integer lenpre,lenpot,ilen,lentmp
2037 character*3 out1file_text,ucase
2040 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2041 call getenv_loc("PREFIX",prefix)
2043 call getenv_loc("POT",pot)
2044 call getenv_loc("DIRTMP",tmpdir)
2045 call getenv_loc("CURDIR",curdir)
2046 call getenv_loc("OUT1FILE",out1file_text)
2047 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2048 out1file_text=ucase(out1file_text)
2049 if (out1file_text(1:1).eq."Y") then
2052 out1file=fg_rank.gt.0
2057 if (lentmp.gt.0) then
2058 write (*,'(80(1h!))')
2059 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2060 write (*,'(80(1h!))')
2061 write (*,*)"All output files will be on node /tmp directory."
2063 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2064 if (me.eq.king) then
2065 write (*,*) "The master node is ",nodename
2066 else if (fg_rank.eq.0) then
2067 write (*,*) "I am the CG slave node ",nodename
2069 write (*,*) "I am the FG slave node ",nodename
2072 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2073 lenpre = lentmp+lenpre+1
2075 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2076 C Get the names and open the input files
2077 #if defined(WINIFL) || defined(WINPGI)
2078 open(1,file=pref_orig(:ilen(pref_orig))//
2079 & '.inp',status='old',readonly,shared)
2080 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2081 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2082 C Get parameter filenames and open the parameter files.
2083 call getenv_loc('BONDPAR',bondname)
2084 open (ibond,file=bondname,status='old',readonly,shared)
2085 call getenv_loc('THETPAR',thetname)
2086 open (ithep,file=thetname,status='old',readonly,shared)
2088 call getenv_loc('THETPARPDB',thetname_pdb)
2089 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2091 call getenv_loc('ROTPAR',rotname)
2092 open (irotam,file=rotname,status='old',readonly,shared)
2094 call getenv_loc('ROTPARPDB',rotname_pdb)
2095 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2097 call getenv_loc('TORPAR',torname)
2098 open (itorp,file=torname,status='old',readonly,shared)
2099 call getenv_loc('TORDPAR',tordname)
2100 open (itordp,file=tordname,status='old',readonly,shared)
2101 call getenv_loc('FOURIER',fouriername)
2102 open (ifourier,file=fouriername,status='old',readonly,shared)
2103 call getenv_loc('ELEPAR',elename)
2104 open (ielep,file=elename,status='old',readonly,shared)
2105 call getenv_loc('SIDEPAR',sidename)
2106 open (isidep,file=sidename,status='old',readonly,shared)
2107 #elif (defined CRAY) || (defined AIX)
2108 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2110 c print *,"Processor",myrank," opened file 1"
2111 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2112 c print *,"Processor",myrank," opened file 9"
2113 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2114 C Get parameter filenames and open the parameter files.
2115 call getenv_loc('BONDPAR',bondname)
2116 open (ibond,file=bondname,status='old',action='read')
2117 c print *,"Processor",myrank," opened file IBOND"
2118 call getenv_loc('THETPAR',thetname)
2119 open (ithep,file=thetname,status='old',action='read')
2120 c print *,"Processor",myrank," opened file ITHEP"
2122 call getenv_loc('THETPARPDB',thetname_pdb)
2123 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2125 call getenv_loc('ROTPAR',rotname)
2126 open (irotam,file=rotname,status='old',action='read')
2127 c print *,"Processor",myrank," opened file IROTAM"
2129 call getenv_loc('ROTPARPDB',rotname_pdb)
2130 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2132 call getenv_loc('TORPAR',torname)
2133 open (itorp,file=torname,status='old',action='read')
2134 c print *,"Processor",myrank," opened file ITORP"
2135 call getenv_loc('TORDPAR',tordname)
2136 open (itordp,file=tordname,status='old',action='read')
2137 c print *,"Processor",myrank," opened file ITORDP"
2138 call getenv_loc('SCCORPAR',sccorname)
2139 open (isccor,file=sccorname,status='old',action='read')
2140 c print *,"Processor",myrank," opened file ISCCOR"
2141 call getenv_loc('FOURIER',fouriername)
2142 open (ifourier,file=fouriername,status='old',action='read')
2143 c print *,"Processor",myrank," opened file IFOURIER"
2144 call getenv_loc('ELEPAR',elename)
2145 open (ielep,file=elename,status='old',action='read')
2146 c print *,"Processor",myrank," opened file IELEP"
2147 call getenv_loc('SIDEPAR',sidename)
2148 open (isidep,file=sidename,status='old',action='read')
2149 c print *,"Processor",myrank," opened file ISIDEP"
2150 c print *,"Processor",myrank," opened parameter files"
2152 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2153 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2154 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2155 C Get parameter filenames and open the parameter files.
2156 call getenv_loc('BONDPAR',bondname)
2157 open (ibond,file=bondname,status='old')
2158 call getenv_loc('THETPAR',thetname)
2159 open (ithep,file=thetname,status='old')
2161 call getenv_loc('THETPARPDB',thetname_pdb)
2162 open (ithep_pdb,file=thetname_pdb,status='old')
2164 call getenv_loc('ROTPAR',rotname)
2165 open (irotam,file=rotname,status='old')
2167 call getenv_loc('ROTPARPDB',rotname_pdb)
2168 open (irotam_pdb,file=rotname_pdb,status='old')
2170 call getenv_loc('TORPAR',torname)
2171 open (itorp,file=torname,status='old')
2172 call getenv_loc('TORDPAR',tordname)
2173 open (itordp,file=tordname,status='old')
2174 call getenv_loc('SCCORPAR',sccorname)
2175 open (isccor,file=sccorname,status='old')
2176 call getenv_loc('FOURIER',fouriername)
2177 open (ifourier,file=fouriername,status='old')
2178 call getenv_loc('ELEPAR',elename)
2179 open (ielep,file=elename,status='old')
2180 call getenv_loc('SIDEPAR',sidename)
2181 open (isidep,file=sidename,status='old')
2183 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2185 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2186 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2187 C Get parameter filenames and open the parameter files.
2188 call getenv_loc('BONDPAR',bondname)
2189 open (ibond,file=bondname,status='old',action='read')
2190 call getenv_loc('THETPAR',thetname)
2191 open (ithep,file=thetname,status='old',action='read')
2193 call getenv_loc('THETPARPDB',thetname_pdb)
2194 print *,"thetname_pdb ",thetname_pdb
2195 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2196 print *,ithep_pdb," opened"
2198 call getenv_loc('ROTPAR',rotname)
2199 open (irotam,file=rotname,status='old',action='read')
2201 call getenv_loc('ROTPARPDB',rotname_pdb)
2202 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2204 call getenv_loc('TORPAR',torname)
2205 open (itorp,file=torname,status='old',action='read')
2206 call getenv_loc('TORDPAR',tordname)
2207 open (itordp,file=tordname,status='old',action='read')
2208 call getenv_loc('SCCORPAR',sccorname)
2209 open (isccor,file=sccorname,status='old',action='read')
2210 call getenv_loc('FOURIER',fouriername)
2211 open (ifourier,file=fouriername,status='old',action='read')
2212 call getenv_loc('ELEPAR',elename)
2213 open (ielep,file=elename,status='old',action='read')
2214 call getenv_loc('SIDEPAR',sidename)
2215 open (isidep,file=sidename,status='old',action='read')
2219 C 8/9/01 In the newest version SCp interaction constants are read from a file
2220 C Use -DOLDSCP to use hard-coded constants instead.
2222 call getenv_loc('SCPPAR',scpname)
2223 #if defined(WINIFL) || defined(WINPGI)
2224 open (iscpp,file=scpname,status='old',readonly,shared)
2225 #elif (defined CRAY) || (defined AIX)
2226 open (iscpp,file=scpname,status='old',action='read')
2228 open (iscpp,file=scpname,status='old')
2230 open (iscpp,file=scpname,status='old',action='read')
2233 call getenv_loc('PATTERN',patname)
2234 #if defined(WINIFL) || defined(WINPGI)
2235 open (icbase,file=patname,status='old',readonly,shared)
2236 #elif (defined CRAY) || (defined AIX)
2237 open (icbase,file=patname,status='old',action='read')
2239 open (icbase,file=patname,status='old')
2241 open (icbase,file=patname,status='old',action='read')
2244 C Open output file only for CG processes
2245 c print *,"Processor",myrank," fg_rank",fg_rank
2246 if (fg_rank.eq.0) then
2248 if (nodes.eq.1) then
2251 npos = dlog10(dfloat(nodes-1))+1
2253 if (npos.lt.3) npos=3
2254 write (liczba,'(i1)') npos
2255 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2257 write (liczba,form) me
2258 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2259 & liczba(:ilen(liczba))
2260 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2262 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2264 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2265 & liczba(:ilen(liczba))//'.mol2'
2266 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2267 & liczba(:ilen(liczba))//'.stat'
2269 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2270 & //liczba(:ilen(liczba))//'.stat')
2271 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2274 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2275 & liczba(:ilen(liczba))//'.const'
2280 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2281 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2282 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2283 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2284 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2286 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2288 rest2name=prefix(:ilen(prefix))//'.rst'
2290 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2293 #if defined(AIX) || defined(PGI)
2294 if (me.eq.king .or. .not. out1file)
2295 & open(iout,file=outname,status='unknown')
2298 if (fg_rank.gt.0) then
2299 write (liczba,'(i3.3)') myrank/nfgtasks
2300 write (ll,'(bz,i3.3)') fg_rank
2301 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2307 open(igeom,file=intname,status='unknown',position='append')
2308 open(ipdb,file=pdbname,status='unknown')
2309 open(imol2,file=mol2name,status='unknown')
2310 open(istat,file=statname,status='unknown',position='append')
2312 c1out open(iout,file=outname,status='unknown')
2315 if (me.eq.king .or. .not.out1file)
2316 & open(iout,file=outname,status='unknown')
2319 if (fg_rank.gt.0) then
2320 print "Processor",fg_rank," opening output file"
2321 write (liczba,'(i3.3)') myrank/nfgtasks
2322 write (ll,'(bz,i3.3)') fg_rank
2323 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2329 open(igeom,file=intname,status='unknown',access='append')
2330 open(ipdb,file=pdbname,status='unknown')
2331 open(imol2,file=mol2name,status='unknown')
2332 open(istat,file=statname,status='unknown',access='append')
2334 c1out open(iout,file=outname,status='unknown')
2337 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2338 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2339 csa csa_history=prefix(:lenpre)//'.CSA.history'
2340 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2341 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2342 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2343 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2344 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2345 csa csa_int=prefix(:lenpre)//'.int'
2346 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2347 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2348 csa csa_in=prefix(:lenpre)//'.CSA.in'
2349 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2352 write (iout,'(80(1h-))')
2353 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2354 write (iout,'(80(1h-))')
2355 write (iout,*) "Input file : ",
2356 & pref_orig(:ilen(pref_orig))//'.inp'
2357 write (iout,*) "Output file : ",
2358 & outname(:ilen(outname))
2360 write (iout,*) "Sidechain potential file : ",
2361 & sidename(:ilen(sidename))
2363 write (iout,*) "SCp potential file : ",
2364 & scpname(:ilen(scpname))
2366 write (iout,*) "Electrostatic potential file : ",
2367 & elename(:ilen(elename))
2368 write (iout,*) "Cumulant coefficient file : ",
2369 & fouriername(:ilen(fouriername))
2370 write (iout,*) "Torsional parameter file : ",
2371 & torname(:ilen(torname))
2372 write (iout,*) "Double torsional parameter file : ",
2373 & tordname(:ilen(tordname))
2374 write (iout,*) "SCCOR parameter file : ",
2375 & sccorname(:ilen(sccorname))
2376 write (iout,*) "Bond & inertia constant file : ",
2377 & bondname(:ilen(bondname))
2378 write (iout,*) "Bending parameter file : ",
2379 & thetname(:ilen(thetname))
2380 write (iout,*) "Rotamer parameter file : ",
2381 & rotname(:ilen(rotname))
2382 write (iout,*) "Threading database : ",
2383 & patname(:ilen(patname))
2385 &write (iout,*)" DIRTMP : ",
2387 write (iout,'(80(1h-))')
2391 c----------------------------------------------------------------------------
2392 subroutine card_concat(card)
2393 implicit real*8 (a-h,o-z)
2394 include 'DIMENSIONS'
2395 include 'COMMON.IOUNITS'
2397 character*80 karta,ucase
2399 read (inp,'(a)') karta
2402 do while (karta(80:80).eq.'&')
2403 card=card(:ilen(card)+1)//karta(:79)
2404 read (inp,'(a)') karta
2407 card=card(:ilen(card)+1)//karta
2410 c----------------------------------------------------------------------------------
2412 implicit real*8 (a-h,o-z)
2413 include 'DIMENSIONS'
2414 include 'COMMON.CHAIN'
2415 include 'COMMON.IOUNITS'
2417 include 'COMMON.CONTROL'
2418 open(irest2,file=rest2name,status='unknown')
2419 read(irest2,*) totT,EK,potE,totE,t_bath
2421 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2424 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2426 if(usampl.or.homol_nset.gt.1) then
2427 read (irest2,*) iset
2432 c---------------------------------------------------------------------------------
2433 subroutine read_fragments
2434 implicit real*8 (a-h,o-z)
2435 include 'DIMENSIONS'
2439 include 'COMMON.SETUP'
2440 include 'COMMON.CHAIN'
2441 include 'COMMON.IOUNITS'
2443 include 'COMMON.CONTROL'
2444 read(inp,*) nset,nfrag,npair,nfrag_back
2445 if(me.eq.king.or..not.out1file)
2446 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2447 & " nfrag_back",nfrag_back
2449 read(inp,*) mset(iset)
2451 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2453 if(me.eq.king.or..not.out1file)
2454 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2455 & ifrag(2,i,iset), qinfrag(i,iset)
2458 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2460 if(me.eq.king.or..not.out1file)
2461 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2462 & ipair(2,i,iset), qinpair(i,iset)
2465 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2466 & wfrag_back(3,i,iset),
2467 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2468 if(me.eq.king.or..not.out1file)
2469 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2470 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2475 c-------------------------------------------------------------------------------
2476 subroutine read_dist_constr
2477 implicit real*8 (a-h,o-z)
2478 include 'DIMENSIONS'
2482 include 'COMMON.SETUP'
2483 include 'COMMON.CONTROL'
2484 include 'COMMON.CHAIN'
2485 include 'COMMON.IOUNITS'
2486 include 'COMMON.SBRIDGE'
2487 integer ifrag_(2,100),ipair_(2,100)
2488 double precision wfrag_(100),wpair_(100)
2489 character*500 controlcard
2490 c write (iout,*) "Calling read_dist_constr"
2491 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2493 call card_concat(controlcard)
2494 call readi(controlcard,"NFRAG",nfrag_,0)
2495 call readi(controlcard,"NPAIR",npair_,0)
2496 call readi(controlcard,"NDIST",ndist_,0)
2497 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2498 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2499 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2500 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2501 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2502 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2503 c write (iout,*) "IFRAG"
2505 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2507 c write (iout,*) "IPAIR"
2509 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2511 if (.not.refstr .and. nfrag.gt.0) then
2513 & "ERROR: no reference structure to compute distance restraints"
2515 & "Restraints must be specified explicitly (NDIST=number)"
2518 if (nfrag.lt.2 .and. npair.gt.0) then
2519 write (iout,*) "ERROR: Less than 2 fragments specified",
2520 & " but distance restraints between pairs requested"
2525 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2526 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2527 & ifrag_(2,i)=nstart_sup+nsup-1
2528 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2530 if (wfrag_(i).gt.0.0d0) then
2531 do j=ifrag_(1,i),ifrag_(2,i)-1
2532 do k=j+1,ifrag_(2,i)
2533 c write (iout,*) "j",j," k",k
2535 if (constr_dist.eq.1) then
2540 forcon(nhpb)=wfrag_(i)
2541 else if (constr_dist.eq.2) then
2542 if (ddjk.le.dist_cut) then
2547 forcon(nhpb)=wfrag_(i)
2554 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2557 if (.not.out1file .or. me.eq.king)
2558 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2559 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2561 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2562 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2569 if (wpair_(i).gt.0.0d0) then
2577 do j=ifrag_(1,ii),ifrag_(2,ii)
2578 do k=ifrag_(1,jj),ifrag_(2,jj)
2582 forcon(nhpb)=wpair_(i)
2583 dhpb(nhpb)=dist(j,k)
2585 if (.not.out1file .or. me.eq.king)
2586 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2587 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2589 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2590 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2597 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2598 & ibecarb(i),forcon(nhpb+1)
2599 if (forcon(nhpb+1).gt.0.0d0) then
2601 if (ibecarb(i).gt.0) then
2602 ihpb(i)=ihpb(i)+nres
2603 jhpb(i)=jhpb(i)+nres
2605 if (dhpb(nhpb).eq.0.0d0)
2606 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2610 if (.not.out1file .or. me.eq.king) then
2613 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2614 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2622 c-------------------------------------------------------------------------------
2624 subroutine read_constr_homology
2626 include 'DIMENSIONS'
2630 include 'COMMON.SETUP'
2631 include 'COMMON.CONTROL'
2632 include 'COMMON.CHAIN'
2633 include 'COMMON.IOUNITS'
2635 include 'COMMON.GEO'
2636 include 'COMMON.INTERACT'
2637 double precision odl_temp,sigma_odl_temp
2638 common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2639 & sigma_odl_temp(maxres,maxres,max_template)
2641 character*24 model_ki_dist, model_ki_angle
2642 character*500 controlcard
2643 character*3200 controlcard1
2644 integer ki, i, j, k, l
2645 logical lprn /.true./
2647 call card_concat(controlcard)
2648 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2649 if (homol_nset.gt.1)then
2650 call card_concat(controlcard)
2651 read(controlcard,*) (waga_dist(i),i=1,homol_nset)
2652 call card_concat(controlcard)
2653 read(controlcard,*) (waga_angle(i),i=1,homol_nset)
2654 write(iout,*) "iset distance_weight angle_weight"
2656 write(iout,*) i,waga_dist(i),waga_angle(i)
2660 call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
2661 call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
2664 write (iout,*) "nnt",nnt," nct",nct
2670 do ki=1,constr_homology
2671 sigma_odl_temp(i,j,ki)=0.0d0
2672 odl_temp(i,j,ki)=0.0d0
2677 do ki=1,constr_homology
2679 sigma_dih(ki,i)=0.0d0
2682 do ki=1,constr_homology
2683 write(kic2,'(i2)') ki
2684 if (ki.le.9) kic2="0"//kic2(2:2)
2686 model_ki_dist="model"//kic2//".dist"
2687 model_ki_angle="model"//kic2//".angle"
2688 open (ientin,file=model_ki_dist,status='old')
2689 do irec=1,maxdim !petla do czytania wiezow na odleglosc
2690 read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
2691 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
2692 odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
2693 sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
2694 & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
2698 open (ientin,file=model_ki_angle,status='old')
2699 do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
2700 read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
2701 & sigma_dih(ki,i+nnt-1)
2702 if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
2703 sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
2709 write (iout,*) "nnt",nnt," nct",nct
2713 c write (iout,*) "i",i," j",j," constr_homology",constr_homology
2714 do while (ki.le.constr_homology .and.
2715 & sigma_odl_temp(i,j,ki).le.0.0d0)
2716 c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
2719 c write (iout,*) "ki",ki
2720 if (ki.gt.constr_homology) cycle
2724 do ki=1,constr_homology
2725 odl(ki,ii)=odl_temp(i,j,ki)
2726 sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
2731 if (constr_homology.gt.0) call homology_partition
2733 if (.not.lprn) return
2734 write (iout,*) "Distance restraints from templates"
2736 write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
2737 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
2739 write (iout,*) "Dihedral angle restraints from templates"
2741 write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2742 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2744 c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
2745 c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
2750 c----------------------------------------------------------------------
2753 subroutine flush(iu)
2758 subroutine flush(iu)
2764 c------------------------------------------------------------------------------
2765 subroutine copy_to_tmp(source)
2766 include "DIMENSIONS"
2767 include "COMMON.IOUNITS"
2768 character*(*) source
2769 character* 256 tmpfile
2773 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2774 inquire(file=tmpfile,exist=ex)
2776 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2777 & " to temporary directory..."
2778 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2779 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2783 c------------------------------------------------------------------------------
2784 subroutine move_from_tmp(source)
2785 include "DIMENSIONS"
2786 include "COMMON.IOUNITS"
2787 character*(*) source
2790 write (*,*) "Moving ",source(:ilen(source)),
2791 & " from temporary directory to working directory"
2792 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2793 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2796 c------------------------------------------------------------------------------
2797 subroutine random_init(seed)
2799 C Initialize random number generator
2801 implicit real*8 (a-h,o-z)
2802 include 'DIMENSIONS'
2808 logical OKRandom, prng_restart
2810 integer iseed_array(4)
2812 include 'COMMON.IOUNITS'
2813 include 'COMMON.TIME1'
2814 include 'COMMON.THREAD'
2815 include 'COMMON.SBRIDGE'
2816 include 'COMMON.CONTROL'
2817 include 'COMMON.MCM'
2818 include 'COMMON.MAP'
2819 include 'COMMON.HEADER'
2820 csa include 'COMMON.CSA'
2821 include 'COMMON.CHAIN'
2822 include 'COMMON.MUCA'
2824 include 'COMMON.FFIELD'
2825 include 'COMMON.SETUP'
2826 iseed=-dint(dabs(seed))
2827 if (iseed.eq.0) then
2828 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2829 & 'Random seed undefined. The program will stop.'
2830 write (*,'(/80(1h*)/20x,a/80(1h*))')
2831 & 'Random seed undefined. The program will stop.'
2833 call mpi_finalize(mpi_comm_world,ierr)
2835 stop 'Bad random seed.'
2838 if (fg_rank.eq.0) then
2842 if(me.eq.king .or. .not. out1file)
2843 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2844 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2845 OKRandom = prng_restart(me,iseedi8)
2848 tmp=65536.0d0**(4-i)
2849 iseed_array(i) = dint(seed/tmp)
2850 seed=seed-iseed_array(i)*tmp
2852 if(me.eq.king .or. .not. out1file)
2853 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2854 & (iseed_array(i),i=1,4)
2855 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2856 & (iseed_array(i),i=1,4)
2857 OKRandom = prng_restart(me,iseed_array)
2860 r1=ran_number(0.0D0,1.0D0)
2861 if(me.eq.king .or. .not. out1file)
2862 & write (iout,*) 'ran_num',r1
2863 if (r1.lt.0.0d0) OKRandom=.false.
2865 if (.not.OKRandom) then
2866 write (iout,*) 'PRNG IS NOT WORKING!!!'
2867 print *,'PRNG IS NOT WORKING!!!'
2870 call mpi_abort(mpi_comm_world,error_msg,ierr)
2873 write (iout,*) 'too many processors for parallel prng'
2874 write (*,*) 'too many processors for parallel prng'
2882 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)