2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read job setup parameters
15 C Read force-field parameters except weights
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 c Read energy-term weights and disulfide parameters
39 C Read molecule information, molecule geometry, energy-term weights, and
40 C restraints if requested
42 C Print restraint information
44 if (.not. out1file .or. me.eq.king) then
47 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
48 & "The following",nhpb-nss,
49 & " distance restraints have been imposed:",
50 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
53 write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
54 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
59 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
60 & "The following",npeak,
61 & " NMR peak restraints have been imposed:",
62 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
63 & " score"," type"," ipeak"
65 do j=ipeak(1,i),ipeak(2,i)
66 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
67 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
68 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
71 write (iout,*) "The ipeak array"
73 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
79 c print *,"Processor",myrank," leaves READRTNS"
82 C-------------------------------------------------------------------------------
83 subroutine read_control
87 implicit real*8 (a-h,o-z)
91 logical OKRandom, prng_restart
94 include 'COMMON.IOUNITS'
95 include 'COMMON.TIME1'
96 include 'COMMON.THREAD'
97 include 'COMMON.SBRIDGE'
98 include 'COMMON.CONTROL'
101 include 'COMMON.HEADER'
103 include 'COMMON.CHAIN'
104 include 'COMMON.MUCA'
106 include 'COMMON.FFIELD'
107 include 'COMMON.INTERACT'
108 include 'COMMON.SETUP'
109 include 'COMMON.SPLITELE'
110 include 'COMMON.SHIELD'
112 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
113 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
115 character*320 controlcard
120 read (INP,'(a)') titel
121 call card_concat(controlcard)
122 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
123 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
124 call reada(controlcard,'SEED',seed,0.0D0)
125 call random_init(seed)
126 C Set up the time limit (caution! The time must be input in minutes!)
127 read_cart=index(controlcard,'READ_CART').gt.0
128 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
129 C this variable with_theta_constr is the variable which allow to read and execute the
130 C constrains on theta angles WITH_THETA_CONSTR is the keyword
131 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
132 write (iout,*) "with_theta_constr ",with_theta_constr
133 call readi(controlcard,'NSAXS',nsaxs,0)
134 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
135 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
136 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
137 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
138 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
139 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
140 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
141 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
142 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
143 call readi(controlcard,'SYM',symetr,1)
144 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
145 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
146 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
147 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
148 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
149 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
150 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
151 call reada(controlcard,'DRMS',drms,0.1D0)
152 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
153 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
154 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
155 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
156 write (iout,'(a,f10.1)')'DRMS = ',drms
157 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
158 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
160 call readi(controlcard,'NZ_START',nz_start,0)
161 call readi(controlcard,'NZ_END',nz_end,0)
162 c call readi(controlcard,'IZ_SC',iz_sc,0)
164 safety = 60.0d0*safety
167 call reada(controlcard,"T_BATH",t_bath,300.0d0)
168 minim=(index(controlcard,'MINIMIZE').gt.0)
169 dccart=(index(controlcard,'CART').gt.0)
170 overlapsc=(index(controlcard,'OVERLAP').gt.0)
171 overlapsc=.not.overlapsc
172 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
173 searchsc=.not.searchsc
174 sideadd=(index(controlcard,'SIDEADD').gt.0)
175 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
176 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
177 outpdb=(index(controlcard,'PDBOUT').gt.0)
178 outmol2=(index(controlcard,'MOL2OUT').gt.0)
179 pdbref=(index(controlcard,'PDBREF').gt.0)
180 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
181 indpdb=index(controlcard,'PDBSTART')
182 extconf=(index(controlcard,'EXTCONF').gt.0)
183 AFMlog=(index(controlcard,'AFM'))
184 selfguide=(index(controlcard,'SELFGUIDE'))
185 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
186 call readi(controlcard,'TUBEMOD',tubelog,0)
187 c write (iout,*) TUBElog,"TUBEMODE"
188 call readi(controlcard,'IPRINT',iprint,0)
189 C SHIELD keyword sets if the shielding effect of side-chains is used
190 C 0 denots no shielding is used all peptide are equally despite the
191 C solvent accesible area
192 C 1 the newly introduced function
193 C 2 reseved for further possible developement
194 call readi(controlcard,'SHIELD',shield_mode,0)
195 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
196 write(iout,*) "shield_mode",shield_mode
198 call readi(controlcard,'TORMODE',tor_mode,0)
199 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
200 write(iout,*) "torsional and valence angle mode",tor_mode
201 call readi(controlcard,'MAXGEN',maxgen,10000)
202 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
203 call readi(controlcard,"KDIAG",kdiag,0)
204 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
205 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
206 & write (iout,*) "RESCALE_MODE",rescale_mode
207 split_ene=index(controlcard,'SPLIT_ENE').gt.0
208 if (index(controlcard,'REGULAR').gt.0.0D0) then
209 call reada(controlcard,'WEIDIS',weidis,0.1D0)
214 if (index(controlcard,"CASC").gt.0) then
216 c else if (index(controlcard,"CAONLY").gt.0) then
218 else if (index(controlcard,"SCONLY").gt.0) then
222 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
223 c & " specifying REFSTR, PDBREF or REGULAR."
227 if (index(controlcard,'CHECKGRAD').gt.0) then
229 if (index(controlcard,' CART').gt.0) then
231 elseif (index(controlcard,'CARINT').gt.0) then
236 call reada(controlcard,'DELTA',aincr,1.0d-7)
237 c write (iout,*) "icheckgrad",icheckgrad
238 elseif (index(controlcard,'THREAD').gt.0) then
240 call readi(controlcard,'THREAD',nthread,0)
241 if (nthread.gt.0) then
242 call reada(controlcard,'WEIDIS',weidis,0.1D0)
245 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
246 stop 'Error termination in Read_Control.'
248 else if (index(controlcard,'MCMA').gt.0) then
250 else if (index(controlcard,'MCEE').gt.0) then
252 else if (index(controlcard,'MULTCONF').gt.0) then
254 else if (index(controlcard,'MAP').gt.0) then
256 call readi(controlcard,'MAP',nmap,0)
257 else if (index(controlcard,'CSA').gt.0) then
259 crc else if (index(controlcard,'ZSCORE').gt.0) then
261 crc ZSCORE is rm from UNRES, modecalc=9 is available
264 cfcm else if (index(controlcard,'MCMF').gt.0) then
266 else if (index(controlcard,'SOFTREG').gt.0) then
268 else if (index(controlcard,'CHECK_BOND').gt.0) then
270 else if (index(controlcard,'TEST').gt.0) then
272 else if (index(controlcard,'MD').gt.0) then
274 else if (index(controlcard,'RE ').gt.0) then
278 lmuca=index(controlcard,'MUCA').gt.0
279 call readi(controlcard,'MUCADYN',mucadyn,0)
280 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
281 if (lmuca .and. (me.eq.king .or. .not.out1file ))
283 write (iout,*) 'MUCADYN=',mucadyn
284 write (iout,*) 'MUCASMOOTH=',muca_smooth
287 iscode=index(controlcard,'ONE_LETTER')
288 indphi=index(controlcard,'PHI')
289 indback=index(controlcard,'BACK')
290 iranconf=index(controlcard,'RAND_CONF')
291 i2ndstr=index(controlcard,'USE_SEC_PRED')
292 gradout=index(controlcard,'GRADOUT').gt.0
293 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
294 C DISTCHAINMAX become obsolete for periodic boundry condition
295 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
296 C Reading the dimensions of box in x,y,z coordinates
297 call reada(controlcard,'BOXX',boxxsize,100.0d0)
298 call reada(controlcard,'BOXY',boxysize,100.0d0)
299 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
300 c Cutoff range for interactions
301 call reada(controlcard,"R_CUT",r_cut,15.0d0)
302 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
303 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
304 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
305 if (lipthick.gt.0.0d0) then
306 bordliptop=(boxzsize+lipthick)/2.0
307 bordlipbot=bordliptop-lipthick
309 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
310 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
311 buflipbot=bordlipbot+lipbufthick
312 bufliptop=bordliptop-lipbufthick
313 if ((lipbufthick*2.0d0).gt.lipthick)
314 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
316 write(iout,*) "bordliptop=",bordliptop
317 write(iout,*) "bordlipbot=",bordlipbot
318 write(iout,*) "bufliptop=",bufliptop
319 write(iout,*) "buflipbot=",buflipbot
320 write (iout,*) "SHIELD MODE",shield_mode
321 if (TUBElog.gt.0) then
322 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
323 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
324 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
325 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
326 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
327 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
328 buftubebot=bordtubebot+tubebufthick
329 buftubetop=bordtubetop-tubebufthick
331 c if (shield_mode.gt.0) then
333 C VSolvSphere the volume of solving sphere
335 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
336 C there will be no distinction between proline peptide group and normal peptide
337 C group in case of shielding parameters
338 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
339 c VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
340 c VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
341 c write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
343 C long axis of side chain
345 c long_r_sidechain(i)=vbldsc0(1,i)
346 c short_r_sidechain(i)=sigma0(i)
350 if (me.eq.king .or. .not.out1file )
351 & write (iout,*) "DISTCHAINMAX",distchainmax
353 if(me.eq.king.or..not.out1file)
354 & write (iout,'(2a)') diagmeth(kdiag),
355 & ' routine used to diagonalize matrices.'
358 c--------------------------------------------------------------------------
359 subroutine read_REMDpar
363 implicit real*8 (a-h,o-z)
365 include 'COMMON.IOUNITS'
366 include 'COMMON.TIME1'
369 include 'COMMON.LANGEVIN'
371 include 'COMMON.LANGEVIN.lang0'
373 include 'COMMON.INTERACT'
374 include 'COMMON.NAMES'
376 include 'COMMON.REMD'
377 include 'COMMON.CONTROL'
378 include 'COMMON.SETUP'
380 character*320 controlcard
381 character*3200 controlcard1
382 integer iremd_m_total
384 if(me.eq.king.or..not.out1file)
385 & write (iout,*) "REMD setup"
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"
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
413 if (index(controlcard,'TLIST').gt.0) then
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)
421 if (index(controlcard,'MLIST').gt.0) then
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)
429 iremd_m_total=iremd_m_total+remd_m(i)
431 write (iout,*) 'Total number of replicas ',iremd_m_total
436 & "Adaptive (PMF-biased) umbrella sampling will be run"
439 if(me.eq.king.or..not.out1file)
440 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
443 c--------------------------------------------------------------------------
444 subroutine read_MDpar
448 implicit real*8 (a-h,o-z)
450 include 'COMMON.IOUNITS'
451 include 'COMMON.TIME1'
454 include 'COMMON.LANGEVIN'
456 include 'COMMON.LANGEVIN.lang0'
458 include 'COMMON.INTERACT'
459 include 'COMMON.NAMES'
461 include 'COMMON.SETUP'
462 include 'COMMON.CONTROL'
463 include 'COMMON.SPLITELE'
464 include 'COMMON.FFIELD'
466 character*320 controlcard
468 call card_concat(controlcard)
469 call readi(controlcard,"NSTEP",n_timestep,1000000)
470 call readi(controlcard,"NTWE",ntwe,100)
471 call readi(controlcard,"NTWX",ntwx,1000)
472 call reada(controlcard,"DT",d_time,1.0d-1)
473 call reada(controlcard,"DVMAX",dvmax,2.0d1)
474 call reada(controlcard,"DAMAX",damax,1.0d1)
475 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
476 call readi(controlcard,"LANG",lang,0)
477 RESPA = index(controlcard,"RESPA") .gt. 0
478 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
479 ntime_split0=ntime_split
480 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
481 ntime_split0=ntime_split
482 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
483 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
484 rest = index(controlcard,"REST").gt.0
485 tbf = index(controlcard,"TBF").gt.0
486 usampl = index(controlcard,"USAMPL").gt.0
487 scale_umb = index(controlcard,"SCALE_UMB").gt.0
488 adaptive = index(controlcard,"ADAPTIVE").gt.0
489 mdpdb = index(controlcard,"MDPDB").gt.0
490 call reada(controlcard,"T_BATH",t_bath,300.0d0)
491 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
492 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
493 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
494 if (count_reset_moment.eq.0) count_reset_moment=1000000000
495 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
496 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
497 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
498 if (count_reset_vel.eq.0) count_reset_vel=1000000000
499 large = index(controlcard,"LARGE").gt.0
500 print_compon = index(controlcard,"PRINT_COMPON").gt.0
501 rattle = index(controlcard,"RATTLE").gt.0
502 preminim = index(controlcard,"PREMINIM").gt.0
503 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
507 write (iout,*) "PREMINIM ",preminim
509 if (index(controlcard,'CART').gt.0) dccart=.true.
510 write (iout,*) "dccart ",dccart
511 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
512 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
514 c if performing umbrella sampling, fragments constrained are read from the fragment file
517 write (iout,*) "Umbrella sampling will be run"
518 if (scale_umb.and.adaptive) then
519 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
520 write (iout,*) "Select one of those and re-run the job."
523 if (scale_umb) write (iout,*)
524 &"Umbrella-restraint force constants will be scaled by temperature"
528 if(me.eq.king.or..not.out1file) then
530 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
532 write (iout,'(a)') "The units are:"
533 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
534 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
535 & " acceleration: angstrom/(48.9 fs)**2"
536 write (iout,'(a)') "energy: kcal/mol, temperature: K"
538 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
539 write (iout,'(a60,f10.5,a)')
540 & "Initial time step of numerical integration:",d_time,
542 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
544 write (iout,'(2a,i4,a)')
545 & "A-MTS algorithm used; initial time step for fast-varying",
546 & " short-range forces split into",ntime_split," steps."
547 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
548 & r_cut," lambda",rlamb
550 write (iout,'(2a,f10.5)')
551 & "Maximum acceleration threshold to reduce the time step",
552 & "/increase split number:",damax
553 write (iout,'(2a,f10.5)')
554 & "Maximum predicted energy drift to reduce the timestep",
555 & "/increase split number:",edriftmax
556 write (iout,'(a60,f10.5)')
557 & "Maximum velocity threshold to reduce velocities:",dvmax
558 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
559 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
560 if (rattle) write (iout,'(a60)')
561 & "Rattle algorithm used to constrain the virtual bonds"
565 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
566 call reada(controlcard,"RWAT",rwat,1.4d0)
567 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
568 surfarea=index(controlcard,"SURFAREA").gt.0
569 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
570 if(me.eq.king.or..not.out1file)then
571 write (iout,'(/a,$)') "Langevin dynamics calculation"
574 & " with direct integration of Langevin equations"
575 else if (lang.eq.2) then
576 write (iout,'(a/)') " with TINKER stochasic MD integrator"
577 else if (lang.eq.3) then
578 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
579 else if (lang.eq.4) then
580 write (iout,'(a/)') " in overdamped mode"
582 write (iout,'(//a,i5)')
583 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
586 write (iout,'(a60,f10.5)') "Temperature:",t_bath
587 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
588 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
589 write (iout,'(a60,f10.5)')
590 & "Scaling factor of the friction forces:",scal_fric
591 if (surfarea) write (iout,'(2a,i10,a)')
592 & "Friction coefficients will be scaled by solvent-accessible",
593 & " surface area every",reset_fricmat," steps."
595 c Calculate friction coefficients and bounds of stochastic forces
596 eta=6*pi*cPoise*etawat
597 if(me.eq.king.or..not.out1file)
598 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
600 gamp=scal_fric*(pstok+rwat)*eta
601 stdfp=dsqrt(2*Rb*t_bath/d_time)
603 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
604 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
606 if(me.eq.king.or..not.out1file)then
607 write (iout,'(/2a/)')
608 & "Radii of site types and friction coefficients and std's of",
609 & " stochastic forces of fully exposed sites"
610 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
612 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
613 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
617 if(me.eq.king.or..not.out1file)then
618 write (iout,'(a)') "Berendsen bath calculation"
619 write (iout,'(a60,f10.5)') "Temperature:",t_bath
620 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
622 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
623 & count_reset_moment," steps"
625 & write (iout,'(a,i10,a)')
626 & "Velocities will be reset at random every",count_reset_vel,
630 if(me.eq.king.or..not.out1file)
631 & write (iout,'(a31)') "Microcanonical mode calculation"
633 if(me.eq.king.or..not.out1file)then
634 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
636 write(iout,*) "MD running with constraints."
637 write(iout,*) "Equilibration time ", eq_time, " mtus."
638 write(iout,*) "Constraining ", nfrag," fragments."
639 write(iout,*) "Length of each fragment, weight and q0:"
641 write (iout,*) "Set of restraints #",iset
643 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
644 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
646 write(iout,*) "constraints between ", npair, "fragments."
647 write(iout,*) "constraint pairs, weights and q0:"
649 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
650 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
652 write(iout,*) "angle constraints within ", nfrag_back,
653 & "backbone fragments."
655 write(iout,*) "fragment, weights, q0:"
657 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
658 & ifrag_back(2,i,iset),
659 & wfrag_back(1,i,iset),qin_back(1,i,iset),
660 & wfrag_back(2,i,iset),qin_back(2,i,iset),
661 & wfrag_back(3,i,iset),qin_back(3,i,iset)
664 write(iout,*) "fragment, weights:"
666 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
667 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
668 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
672 iset=mod(kolor,nset)+1
675 if(me.eq.king.or..not.out1file)
676 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
679 c------------------------------------------------------------------------------
682 C Read molecular data.
684 implicit real*8 (a-h,o-z)
690 include 'COMMON.IOUNITS'
693 include 'COMMON.INTERACT'
694 include 'COMMON.LOCAL'
695 include 'COMMON.NAMES'
696 include 'COMMON.CHAIN'
697 include 'COMMON.FFIELD'
698 include 'COMMON.SBRIDGE'
699 include 'COMMON.HEADER'
700 include 'COMMON.CONTROL'
701 include 'COMMON.DBASE'
702 include 'COMMON.THREAD'
703 include 'COMMON.CONTACTS'
704 include 'COMMON.TORCNSTR'
705 include 'COMMON.TIME1'
706 include 'COMMON.BOUNDS'
708 include 'COMMON.SETUP'
709 include 'COMMON.SHIELD'
710 character*4 sequence(maxres)
712 double precision x(maxvar)
713 character*256 pdbfile
714 character*400 weightcard
715 character*80 weightcard_t,ucase
716 dimension itype_pdb(maxres)
717 common /pizda/ itype_pdb
718 logical seq_comp,fail
719 double precision energia(0:n_ene)
720 double precision secprob(3,maxdih_constr)
725 C Read PDB structure if applicable
727 if (indpdb.gt.0 .or. pdbref) then
728 read(inp,'(a)') pdbfile
729 if(me.eq.king.or..not.out1file)
730 & write (iout,'(2a)') 'PDB data will be read from file ',
731 & pdbfile(:ilen(pdbfile))
732 open(ipdbin,file=pdbfile,status='old',err=33)
734 33 write (iout,'(a)') 'Error opening PDB file.'
745 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
746 & (crefjlee(j,i+nres),j=1,3)
749 if(me.eq.king.or..not.out1file)
750 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
751 & ' nstart_sup=',nstart_sup
753 itype_pdb(i)=itype(i)
757 nct=nstart_sup+nsup-1
758 call contact(.false.,ncont_ref,icont_ref,co)
761 if(me.eq.king.or..not.out1file)
762 & write(iout,*)'Adding sidechains'
766 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
769 do while (fail.and.nsi.le.maxsi)
770 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
773 if(fail) write(iout,*)'Adding sidechain failed for res ',
774 & i,' after ',nsi,' trials'
779 if (indpdb.eq.0) then
780 C Read sequence if not taken from the pdb file.
782 c print *,'nres=',nres
783 if (iscode.gt.0) then
784 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
786 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
788 C Convert sequence to numeric code
790 itype(i)=rescode(i,sequence(i),iscode)
792 C Assign initial virtual bond lengths
797 c if (itype(1).eq.ntyp1) then
799 c vbld_inv(2)=vbld_inv(2)*2
801 c if (itype(nres).eq.ntyp1) then
802 c vbld(nres)=vbld(nres)/2
803 c vbld_inv(nres)=vbld_inv(nres)*2
806 c vbld(i+nres)=dsc(iabs(itype(i)))
807 c vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
808 c write (iout,*) "i",i," itype",itype(i),
809 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
813 c print '(20i4)',(itype(i),i=1,nres)
816 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
818 if (itype(i).eq.ntyp1) then
822 else if (iabs(itype(i+1)).ne.20) then
824 else if (iabs(itype(i)).ne.20) then
831 if(me.eq.king.or..not.out1file)then
832 write (iout,*) "ITEL"
834 write (iout,*) i,itype(i),itel(i)
836 print *,'Call Read_Bridge.'
840 cd print *,'NNT=',NNT,' NCT=',NCT
841 call seq2chains(nres,itype,nchain,chain_length,chain_border,
843 write(iout,*) "nres",nres," nchain",nchain
845 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
848 call chain_symmetry(nchain,nres,itype,chain_border,
849 & chain_length,npermchain,tabpermchain)
851 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
854 write(iout,*) "residue permutations"
856 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
858 if (itype(1).eq.ntyp1) nnt=2
859 if (itype(nres).eq.ntyp1) nct=nct-1
861 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
862 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
864 print*, 'init_dfa_vars finished!'
866 print*, 'read_dfa_info finished!'
870 if(me.eq.king.or..not.out1file)
871 & write (iout,'(a,i3)') 'nsup=',nsup
873 if (nsup.le.(nct-nnt+1)) then
874 do i=0,nct-nnt+1-nsup
875 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
881 & 'Error - sequences to be superposed do not match.'
884 do i=0,nsup-(nct-nnt+1)
885 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
887 nstart_sup=nstart_sup+i
893 & 'Error - sequences to be superposed do not match.'
896 if (nsup.eq.0) nsup=nct-nnt
897 if (nstart_sup.eq.0) nstart_sup=nnt
898 if (nstart_seq.eq.0) nstart_seq=nnt
899 if(me.eq.king.or..not.out1file)
900 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
901 & ' nstart_seq=',nstart_seq
904 C 8/13/98 Set limits to generating the dihedral angles
909 read (inp,*) ndih_constr
910 write (iout,*) "ndih_constr",ndih_constr
911 if (ndih_constr.gt.0) then
914 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
917 if(me.eq.king.or..not.out1file)then
920 & 'There are',ndih_constr,' restraints on gamma angles.'
922 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
929 phi0(i)=deg2rad*phi0(i)
930 drange(i)=deg2rad*drange(i)
932 C if(me.eq.king.or..not.out1file)
933 C & write (iout,*) 'FTORS',ftors
936 phibound(1,ii) = phi0(i)-drange(i)
937 phibound(2,ii) = phi0(i)+drange(i)
940 if (me.eq.king .or. .not.out1file) then
942 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
944 write (iout,'(a3,i5,2f10.1)')
945 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
950 else if (ndih_constr.lt.0) then
952 call card_concat(weightcard)
953 call reada(weightcard,"PHIHEL",phihel,50.0D0)
954 call reada(weightcard,"PHIBET",phibet,180.0D0)
955 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
956 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
957 call reada(weightcard,"WDIHC",wdihc,0.591D0)
958 write (iout,*) "Weight of dihedral angle restraints",wdihc
959 read(inp,'(9x,3f7.3)')
960 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
961 write (iout,*) "The secprob array"
963 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
967 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
968 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
969 ndih_constr=ndih_constr+1
970 idih_constr(ndih_constr)=i
973 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
974 sumv=sumv+vpsipred(j,ndih_constr)
977 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
979 phibound(1,ndih_constr)=phihel*deg2rad
980 phibound(2,ndih_constr)=phibet*deg2rad
981 sdihed(1,ndih_constr)=sigmahel*deg2rad
982 sdihed(2,ndih_constr)=sigmabet*deg2rad
986 if(me.eq.king.or..not.out1file)then
989 & 'There are',ndih_constr,
990 & ' bimodal restraints on gamma angles.'
992 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
993 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
994 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
995 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
996 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
997 & (vpsipred(j,i),j=1,3)
1004 C first setting the theta boundaries to 0 to pi
1005 C this mean that there is no energy penalty for any angle occuring this can be applied
1006 C for generate random conformation but is not implemented in this way
1009 C thetabound(2,i)=pi
1011 C begin reading theta constrains this is quartic constrains allowing to
1012 C have smooth second derivative
1013 if (with_theta_constr) then
1014 C with_theta_constr is keyword allowing for occurance of theta constrains
1015 read (inp,*) ntheta_constr
1016 C ntheta_constr is the number of theta constrains
1017 if (ntheta_constr.gt.0) then
1018 C read (inp,*) ftors
1019 read (inp,*) (itheta_constr(i),theta_constr0(i),
1020 & theta_drange(i),for_thet_constr(i),
1021 & i=1,ntheta_constr)
1022 C the above code reads from 1 to ntheta_constr
1023 C itheta_constr(i) residue i for which is theta_constr
1024 C theta_constr0 the global minimum value
1025 C theta_drange is range for which there is no energy penalty
1026 C for_thet_constr is the force constant for quartic energy penalty
1028 if(me.eq.king.or..not.out1file)then
1030 & 'There are',ntheta_constr,' constraints on phi angles.'
1031 do i=1,ntheta_constr
1032 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1034 & for_thet_constr(i)
1037 do i=1,ntheta_constr
1038 theta_constr0(i)=deg2rad*theta_constr0(i)
1039 theta_drange(i)=deg2rad*theta_drange(i)
1041 C if(me.eq.king.or..not.out1file)
1042 C & write (iout,*) 'FTORS',ftors
1043 C do i=1,ntheta_constr
1044 C ii = itheta_constr(i)
1045 C thetabound(1,ii) = phi0(i)-drange(i)
1046 C thetabound(2,ii) = phi0(i)+drange(i)
1048 endif ! ntheta_constr.gt.0
1049 endif! with_theta_constr
1051 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1052 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1053 c--- Zscore rms -------
1054 if (nz_start.eq.0) nz_start=nnt
1055 if (nz_end.eq.0 .and. nsup.gt.0) then
1057 else if (nz_end.eq.0) then
1060 if(me.eq.king.or..not.out1file)then
1061 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1062 write (iout,*) 'IZ_SC=',iz_sc
1064 c----------------------
1067 if (.not.pdbref) then
1068 call read_angles(inp,*38)
1071 38 write (iout,'(a)') 'Error reading reference structure.'
1073 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1074 stop 'Error reading reference structure'
1076 39 call chainbuild_extconf
1078 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1087 call contact(.true.,ncont_ref,icont_ref,co)
1091 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1093 if (constr_dist.gt.0) call read_dist_constr
1094 write (iout,*) "After read_dist_constr nhpb",nhpb
1095 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1097 call NMRpeak_partition
1098 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1100 if(me.eq.king.or..not.out1file)
1101 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1104 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1106 if(me.eq.king.or..not.out1file)
1107 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1108 & icont_ref(1,i),' ',
1109 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1112 write (iout,*) "calling read_saxs_consrtr",nsaxs
1113 if (nsaxs.gt.0) call read_saxs_constr
1115 if (constr_homology.gt.0) then
1116 call read_constr_homology
1117 if (indpdb.gt.0 .or. pdbref) then
1120 c(j,i)=crefjlee(j,i)
1121 cref(j,i)=crefjlee(j,i)
1126 write (iout,*) "sc_loc_geom: Array C"
1128 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1129 & (c(j,i+nres),j=1,3)
1131 write (iout,*) "Array Cref"
1133 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1134 & (cref(j,i+nres),j=1,3)
1137 call int_from_cart1(.false.)
1138 call sc_loc_geom(.false.)
1140 thetaref(i)=theta(i)
1145 dc(j,i)=c(j,i+1)-c(j,i)
1146 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1151 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1152 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1161 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1162 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1163 & modecalc.ne.10) then
1164 C If input structure hasn't been supplied from the PDB file read or generate
1166 if (iranconf.eq.0 .and. .not. extconf) then
1167 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1168 & write (iout,'(a)') 'Initial geometry will be read in.'
1170 read(inp,'(8f10.5)',end=36,err=36)
1171 & ((c(l,k),l=1,3),k=1,nres),
1172 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1173 write (iout,*) "Exit READ_CART"
1174 c write (iout,'(8f10.5)')
1175 c & ((c(l,k),l=1,3),k=1,nres),
1176 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1180 dc(j,i)=c(j,i+1)-c(j,i)
1181 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1185 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1187 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1188 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1193 c dc_norm(j,i+nres)=0.0d0
1197 call int_from_cart1(.true.)
1198 write (iout,*) "Finish INT_TO_CART"
1199 c write (iout,*) "DC"
1201 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1202 c & (dc(j,i+nres),j=1,3)
1208 call read_angles(inp,*36)
1210 call chainbuild_extconf
1213 36 write (iout,'(a)') 'Error reading angle file.'
1215 call mpi_finalize( MPI_COMM_WORLD,IERR )
1217 stop 'Error reading angle file.'
1219 else if (extconf) then
1220 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1221 & write (iout,'(a)') 'Extended chain initial geometry.'
1223 theta(i)=90d0*deg2rad
1226 phi(i)=180d0*deg2rad
1229 alph(i)=110d0*deg2rad
1232 omeg(i)=-120d0*deg2rad
1233 if (itype(i).le.0) omeg(i)=-omeg(i)
1236 call chainbuild_extconf
1238 if(me.eq.king.or..not.out1file)
1239 & write (iout,'(a)') 'Random-generated initial geometry.'
1242 if (me.eq.king .or. fg_rank.eq.0 .and. (
1243 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1247 call gen_rand_conf(itmp,*30)
1249 30 write (iout,*) 'Failed to generate random conformation',
1250 & ', itrial=',itrial
1251 write (*,*) 'Processor:',me,
1252 & ' Failed to generate random conformation',
1262 write (iout,'(a,i3,a)') 'Processor:',me,
1263 & ' error in generating random conformation.'
1264 write (*,'(a,i3,a)') 'Processor:',me,
1265 & ' error in generating random conformation.'
1268 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1273 & ' error in generating random conformation.'
1278 elseif (modecalc.eq.4) then
1279 read (inp,'(a)') intinname
1280 open (intin,file=intinname,status='old',err=333)
1281 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1282 & write (iout,'(a)') 'intinname',intinname
1283 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1285 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1287 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1289 stop 'Error opening angle file.'
1293 C Generate distance constraints, if the PDB structure is to be regularized.
1294 if (nthread.gt.0) then
1295 call read_threadbase
1298 if (me.eq.king .or. .not. out1file)
1300 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1301 write (iout,'(/a,i3,a)')
1302 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1303 write (iout,'(20i4)') (iss(i),i=1,ns)
1305 write(iout,*)"Running with dynamic disulfide-bond formation"
1307 write (iout,'(/a/)') 'Pre-formed links are:'
1313 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1314 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1320 if (ns.gt.0.and.dyn_ss) then
1324 forcon(i-nss)=forcon(i)
1331 dyn_ss_mask(iss(i))=.true.
1336 c write (iout,*) "DC"
1338 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1339 c & (dc(j,i+nres),j=1,3)
1341 if (i2ndstr.gt.0) call secstrp2dihc
1342 c call geom_to_var(nvar,x)
1343 c call etotal(energia(0))
1344 c call enerprint(energia(0))
1345 c call briefout(0,etot)
1347 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1348 cd write (iout,'(a)') 'Variable list:'
1349 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1351 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1352 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1353 & 'Processor',myrank,': end reading molecular data.'
1358 c--------------------------------------------------------------------------
1359 logical function seq_comp(itypea,itypeb,length)
1361 integer length,itypea(length),itypeb(length)
1364 if (itypea(i).ne.itypeb(i)) then
1372 c-----------------------------------------------------------------------------
1373 subroutine read_bridge
1374 C Read information about disulfide bridges.
1375 implicit real*8 (a-h,o-z)
1376 include 'DIMENSIONS'
1380 include 'COMMON.IOUNITS'
1381 include 'COMMON.GEO'
1382 include 'COMMON.VAR'
1383 include 'COMMON.INTERACT'
1384 include 'COMMON.LOCAL'
1385 include 'COMMON.NAMES'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.FFIELD'
1388 include 'COMMON.SBRIDGE'
1389 include 'COMMON.HEADER'
1390 include 'COMMON.CONTROL'
1391 include 'COMMON.DBASE'
1392 include 'COMMON.THREAD'
1393 include 'COMMON.TIME1'
1394 include 'COMMON.SETUP'
1395 C Read bridging residues.
1396 read (inp,*) ns,(iss(i),i=1,ns)
1398 if(me.eq.king.or..not.out1file)
1399 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1400 C Check whether the specified bridging residues are cystines.
1402 if (itype(iss(i)).ne.1) then
1403 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1404 & 'Do you REALLY think that the residue ',
1405 & restyp(itype(iss(i))),i,
1406 & ' can form a disulfide bridge?!!!'
1407 write (*,'(2a,i3,a)')
1408 & 'Do you REALLY think that the residue ',
1409 & restyp(itype(iss(i))),i,
1410 & ' can form a disulfide bridge?!!!'
1412 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1417 C Read preformed bridges.
1419 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1421 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1424 C Check if the residues involved in bridges are in the specified list of
1425 C bridging residues.
1428 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1429 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1430 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1431 & ' contains residues present in other pairs.'
1432 write (*,'(a,i3,a)') 'Disulfide pair',i,
1433 & ' contains residues present in other pairs.'
1435 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1441 if (ihpb(i).eq.iss(j)) goto 10
1443 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1446 if (jhpb(i).eq.iss(j)) goto 20
1448 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1454 ihpb(i)=ihpb(i)+nres
1455 jhpb(i)=jhpb(i)+nres
1461 c----------------------------------------------------------------------------
1462 subroutine read_x(kanal,*)
1463 implicit real*8 (a-h,o-z)
1464 include 'DIMENSIONS'
1465 include 'COMMON.GEO'
1466 include 'COMMON.VAR'
1467 include 'COMMON.CHAIN'
1468 include 'COMMON.IOUNITS'
1469 include 'COMMON.CONTROL'
1470 include 'COMMON.LOCAL'
1471 include 'COMMON.INTERACT'
1472 c Read coordinates from input
1474 read(kanal,'(8f10.5)',end=10,err=10)
1475 & ((c(l,k),l=1,3),k=1,nres),
1476 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1479 c(j,2*nres)=c(j,nres)
1481 call int_from_cart1(.false.)
1484 dc(j,i)=c(j,i+1)-c(j,i)
1485 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1489 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1491 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1492 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1500 c----------------------------------------------------------------------------
1501 subroutine read_threadbase
1502 implicit real*8 (a-h,o-z)
1503 include 'DIMENSIONS'
1504 include 'COMMON.IOUNITS'
1505 include 'COMMON.GEO'
1506 include 'COMMON.VAR'
1507 include 'COMMON.INTERACT'
1508 include 'COMMON.LOCAL'
1509 include 'COMMON.NAMES'
1510 include 'COMMON.CHAIN'
1511 include 'COMMON.FFIELD'
1512 include 'COMMON.SBRIDGE'
1513 include 'COMMON.HEADER'
1514 include 'COMMON.CONTROL'
1515 include 'COMMON.DBASE'
1516 include 'COMMON.THREAD'
1517 include 'COMMON.TIME1'
1518 C Read pattern database for threading.
1519 read (icbase,*) nseq
1521 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1522 & nres_base(2,i),nres_base(3,i)
1523 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1525 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1526 c & nres_base(2,i),nres_base(3,i)
1527 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1531 if (weidis.eq.0.0D0) weidis=0.1D0
1540 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1541 write (iout,'(a,i5)') 'nexcl: ',nexcl
1542 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1545 c------------------------------------------------------------------------------
1546 subroutine setup_var
1547 implicit real*8 (a-h,o-z)
1548 include 'DIMENSIONS'
1549 include 'COMMON.IOUNITS'
1550 include 'COMMON.GEO'
1551 include 'COMMON.VAR'
1552 include 'COMMON.INTERACT'
1553 include 'COMMON.LOCAL'
1554 include 'COMMON.NAMES'
1555 include 'COMMON.CHAIN'
1556 include 'COMMON.FFIELD'
1557 include 'COMMON.SBRIDGE'
1558 include 'COMMON.HEADER'
1559 include 'COMMON.CONTROL'
1560 include 'COMMON.DBASE'
1561 include 'COMMON.THREAD'
1562 include 'COMMON.TIME1'
1563 C Set up variable list.
1568 write (iout,*) "SETUP_VAR ialph"
1570 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1572 ialph(i,1)=nvar+nside
1576 if (indphi.gt.0) then
1578 else if (indback.gt.0) then
1583 write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1586 c----------------------------------------------------------------------------
1587 subroutine gen_dist_constr
1588 C Generate CA distance constraints.
1589 implicit real*8 (a-h,o-z)
1590 include 'DIMENSIONS'
1591 include 'COMMON.IOUNITS'
1592 include 'COMMON.GEO'
1593 include 'COMMON.VAR'
1594 include 'COMMON.INTERACT'
1595 include 'COMMON.LOCAL'
1596 include 'COMMON.NAMES'
1597 include 'COMMON.CHAIN'
1598 include 'COMMON.FFIELD'
1599 include 'COMMON.SBRIDGE'
1600 include 'COMMON.HEADER'
1601 include 'COMMON.CONTROL'
1602 include 'COMMON.DBASE'
1603 include 'COMMON.THREAD'
1604 include 'COMMON.TIME1'
1605 dimension itype_pdb(maxres)
1606 common /pizda/ itype_pdb
1608 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1609 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1610 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1612 do i=nstart_sup,nstart_sup+nsup-1
1613 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1614 cd & ' seq_pdb', restyp(itype_pdb(i))
1615 do j=i+2,nstart_sup+nsup-1
1617 ihpb(nhpb)=i+nstart_seq-nstart_sup
1618 jhpb(nhpb)=j+nstart_seq-nstart_sup
1620 dhpb(nhpb)=dist(i,j)
1623 cd write (iout,'(a)') 'Distance constraints:'
1628 cd if (ii.gt.nres) then
1633 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1634 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1635 cd & dhpb(i),forcon(i)
1639 c----------------------------------------------------------------------------
1641 implicit real*8 (a-h,o-z)
1642 include 'DIMENSIONS'
1643 include 'COMMON.MAP'
1644 include 'COMMON.IOUNITS'
1645 character*3 angid(4) /'THE','PHI','ALP','OME'/
1646 character*80 mapcard,ucase
1648 read (inp,'(a)') mapcard
1649 mapcard=ucase(mapcard)
1650 if (index(mapcard,'PHI').gt.0) then
1652 else if (index(mapcard,'THE').gt.0) then
1654 else if (index(mapcard,'ALP').gt.0) then
1656 else if (index(mapcard,'OME').gt.0) then
1659 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1660 stop 'Error - illegal variable spec in MAP card.'
1662 call readi (mapcard,'RES1',res1(imap),0)
1663 call readi (mapcard,'RES2',res2(imap),0)
1664 if (res1(imap).eq.0) then
1665 res1(imap)=res2(imap)
1666 else if (res2(imap).eq.0) then
1667 res2(imap)=res1(imap)
1669 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1671 & 'Error - illegal definition of variable group in MAP.'
1672 stop 'Error - illegal definition of variable group in MAP.'
1674 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1675 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1676 call readi(mapcard,'NSTEP',nstep(imap),0)
1677 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1679 & 'Illegal boundary and/or step size specification in MAP.'
1680 stop 'Illegal boundary and/or step size specification in MAP.'
1685 c----------------------------------------------------------------------------
1687 implicit real*8 (a-h,o-z)
1688 include 'DIMENSIONS'
1689 include 'COMMON.IOUNITS'
1690 include 'COMMON.GEO'
1691 include 'COMMON.CSA'
1692 include 'COMMON.BANK'
1693 include 'COMMON.CONTROL'
1695 character*620 mcmcard
1696 call card_concat(mcmcard)
1698 call readi(mcmcard,'NCONF',nconf,50)
1699 call readi(mcmcard,'NADD',nadd,0)
1700 call readi(mcmcard,'JSTART',jstart,1)
1701 call readi(mcmcard,'JEND',jend,1)
1702 call readi(mcmcard,'NSTMAX',nstmax,500000)
1703 call readi(mcmcard,'N0',n0,1)
1704 call readi(mcmcard,'N1',n1,6)
1705 call readi(mcmcard,'N2',n2,4)
1706 call readi(mcmcard,'N3',n3,0)
1707 call readi(mcmcard,'N4',n4,0)
1708 call readi(mcmcard,'N5',n5,0)
1709 call readi(mcmcard,'N6',n6,10)
1710 call readi(mcmcard,'N7',n7,0)
1711 call readi(mcmcard,'N8',n8,0)
1712 call readi(mcmcard,'N9',n9,0)
1713 call readi(mcmcard,'N14',n14,0)
1714 call readi(mcmcard,'N15',n15,0)
1715 call readi(mcmcard,'N16',n16,0)
1716 call readi(mcmcard,'N17',n17,0)
1717 call readi(mcmcard,'N18',n18,0)
1719 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1721 call readi(mcmcard,'NDIFF',ndiff,2)
1722 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1723 call readi(mcmcard,'IS1',is1,1)
1724 call readi(mcmcard,'IS2',is2,8)
1725 call readi(mcmcard,'NRAN0',nran0,4)
1726 call readi(mcmcard,'NRAN1',nran1,2)
1727 call readi(mcmcard,'IRR',irr,1)
1728 call readi(mcmcard,'NSEED',nseed,20)
1729 call readi(mcmcard,'NTOTAL',ntotal,10000)
1730 call reada(mcmcard,'CUT1',cut1,2.0d0)
1731 call reada(mcmcard,'CUT2',cut2,5.0d0)
1732 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1733 call readi(mcmcard,'ICMAX',icmax,3)
1734 call readi(mcmcard,'IRESTART',irestart,0)
1735 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1738 call reada(mcmcard,'DELE',dele,20.0d0)
1739 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1740 call readi(mcmcard,'IREF',iref,0)
1741 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1742 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1743 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1744 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1745 write (iout,*) "NCONF_IN",nconf_in
1748 c----------------------------------------------------------------------------
1749 cfmc subroutine mcmfread
1750 cfmc implicit real*8 (a-h,o-z)
1751 cfmc include 'DIMENSIONS'
1752 cfmc include 'COMMON.MCMF'
1753 cfmc include 'COMMON.IOUNITS'
1754 cfmc include 'COMMON.GEO'
1755 cfmc character*80 ucase
1756 cfmc character*620 mcmcard
1757 cfmc call card_concat(mcmcard)
1759 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1760 cfmc write(iout,*)'MAXRANT=',maxrant
1761 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1762 cfmc write(iout,*)'MAXFAM=',maxfam
1763 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1764 cfmc write(iout,*)'NNET1=',nnet1
1765 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1766 cfmc write(iout,*)'NNET2=',nnet2
1767 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1768 cfmc write(iout,*)'NNET3=',nnet3
1769 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1770 cfmc write(iout,*)'ILASTT=',ilastt
1771 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1772 cfmc write(iout,*)'MAXSTR=',maxstr
1773 cfmc maxstr_f=maxstr/maxfam
1774 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1775 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1776 cfmc write(iout,*)'NMCMF=',nmcmf
1777 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1778 cfmc write(iout,*)'IFOCUS=',ifocus
1779 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1780 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1781 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1782 cfmc write(iout,*)'INTPRT=',intprt
1783 cfmc call readi(mcmcard,'IPRT',iprt,100)
1784 cfmc write(iout,*)'IPRT=',iprt
1785 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1786 cfmc write(iout,*)'IMAXTR=',imaxtr
1787 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1788 cfmc write(iout,*)'MAXEVEN=',maxeven
1789 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1790 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1791 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1792 cfmc write(iout,*)'INIMIN=',inimin
1793 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1794 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1795 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1796 cfmc write(iout,*)'NTHREAD=',nthread
1797 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1798 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1799 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1800 cfmc write(iout,*)'MAXPERT=',maxpert
1801 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1802 cfmc write(iout,*)'IRMSD=',irmsd
1803 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1804 cfmc write(iout,*)'DENEMIN=',denemin
1805 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1806 cfmc write(iout,*)'RCUT1S=',rcut1s
1807 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1808 cfmc write(iout,*)'RCUT1E=',rcut1e
1809 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1810 cfmc write(iout,*)'RCUT2S=',rcut2s
1811 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1812 cfmc write(iout,*)'RCUT2E=',rcut2e
1813 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1814 cfmc write(iout,*)'DPERT1=',d_pert1
1815 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1816 cfmc write(iout,*)'DPERT1A=',d_pert1a
1817 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1818 cfmc write(iout,*)'DPERT2=',d_pert2
1819 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1820 cfmc write(iout,*)'DPERT2A=',d_pert2a
1821 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1822 cfmc write(iout,*)'DPERT2B=',d_pert2b
1823 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1824 cfmc write(iout,*)'DPERT2C=',d_pert2c
1825 cfmc d_pert1=deg2rad*d_pert1
1826 cfmc d_pert1a=deg2rad*d_pert1a
1827 cfmc d_pert2=deg2rad*d_pert2
1828 cfmc d_pert2a=deg2rad*d_pert2a
1829 cfmc d_pert2b=deg2rad*d_pert2b
1830 cfmc d_pert2c=deg2rad*d_pert2c
1831 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1832 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1833 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1834 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1835 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1836 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1837 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1838 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1839 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1840 cfmc write(iout,*)'RCUTINI=',rcutini
1841 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1842 cfmc write(iout,*)'GRAT=',grat
1843 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1844 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1848 c----------------------------------------------------------------------------
1850 implicit real*8 (a-h,o-z)
1851 include 'DIMENSIONS'
1852 include 'COMMON.MCM'
1853 include 'COMMON.MCE'
1854 include 'COMMON.IOUNITS'
1856 character*320 mcmcard
1857 call card_concat(mcmcard)
1858 call readi(mcmcard,'MAXACC',maxacc,100)
1859 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1860 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1861 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1862 call readi(mcmcard,'MAXREPM',maxrepm,200)
1863 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1864 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1865 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1866 call reada(mcmcard,'E_UP',e_up,5.0D0)
1867 call reada(mcmcard,'DELTE',delte,0.1D0)
1868 call readi(mcmcard,'NSWEEP',nsweep,5)
1869 call readi(mcmcard,'NSTEPH',nsteph,0)
1870 call readi(mcmcard,'NSTEPC',nstepc,0)
1871 call reada(mcmcard,'TMIN',tmin,298.0D0)
1872 call reada(mcmcard,'TMAX',tmax,298.0D0)
1873 call readi(mcmcard,'NWINDOW',nwindow,0)
1874 call readi(mcmcard,'PRINT_MC',print_mc,0)
1875 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1876 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1877 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1878 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1879 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1880 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1881 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1882 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1883 if (nwindow.gt.0) then
1884 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1886 winlen(i)=winend(i)-winstart(i)+1
1889 if (tmax.lt.tmin) tmax=tmin
1890 if (tmax.eq.tmin) then
1894 if (nstepc.gt.0 .and. nsteph.gt.0) then
1895 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1896 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1898 C Probabilities of different move types
1899 sumpro_type(0)=0.0D0
1900 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1901 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1902 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1903 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1904 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1905 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1906 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1908 print *,'i',i,' sumprotype',sumpro_type(i)
1909 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1910 print *,'i',i,' sumprotype',sumpro_type(i)
1914 c----------------------------------------------------------------------------
1915 subroutine read_minim
1916 implicit real*8 (a-h,o-z)
1917 include 'DIMENSIONS'
1918 include 'COMMON.MINIM'
1919 include 'COMMON.IOUNITS'
1921 character*320 minimcard
1922 call card_concat(minimcard)
1923 call readi(minimcard,'MAXMIN',maxmin,2000)
1924 call readi(minimcard,'MAXFUN',maxfun,5000)
1925 call readi(minimcard,'MINMIN',minmin,maxmin)
1926 call readi(minimcard,'MINFUN',minfun,maxmin)
1927 call reada(minimcard,'TOLF',tolf,1.0D-2)
1928 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1929 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1930 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1931 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1932 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1933 & 'Options in energy minimization:'
1934 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1935 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1936 & 'MinMin:',MinMin,' MinFun:',MinFun,
1937 & ' TolF:',TolF,' RTolF:',RTolF
1940 c----------------------------------------------------------------------------
1941 subroutine read_angles(kanal,*)
1942 implicit real*8 (a-h,o-z)
1943 include 'DIMENSIONS'
1944 include 'COMMON.GEO'
1945 include 'COMMON.VAR'
1946 include 'COMMON.CHAIN'
1947 include 'COMMON.IOUNITS'
1948 include 'COMMON.CONTROL'
1949 c Read angles from input
1951 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1952 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1953 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1954 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1957 c 9/7/01 avoid 180 deg valence angle
1958 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1960 theta(i)=deg2rad*theta(i)
1961 phi(i)=deg2rad*phi(i)
1962 alph(i)=deg2rad*alph(i)
1963 omeg(i)=deg2rad*omeg(i)
1968 c----------------------------------------------------------------------------
1969 subroutine reada(rekord,lancuch,wartosc,default)
1971 character*(*) rekord,lancuch
1972 double precision 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 readi(rekord,lancuch,wartosc,default)
1989 character*(*) rekord,lancuch
1990 integer wartosc,default
1993 iread=index(rekord,lancuch)
1994 if (iread.eq.0) then
1998 iread=iread+ilen(lancuch)+1
1999 read (rekord(iread:),*,err=10,end=10) wartosc
2004 c----------------------------------------------------------------------------
2005 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2008 integer 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 multreada(rekord,lancuch,tablica,dim,default)
2026 double precision tablica(dim),default
2027 character*(*) rekord,lancuch
2034 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2035 if (iread.eq.0) return
2036 iread=iread+ilen(lancuch)+1
2037 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2040 c----------------------------------------------------------------------------
2041 subroutine openunits
2042 implicit real*8 (a-h,o-z)
2043 include 'DIMENSIONS'
2046 character*16 form,nodename
2049 include 'COMMON.SETUP'
2050 include 'COMMON.IOUNITS'
2052 include 'COMMON.CONTROL'
2053 integer lenpre,lenpot,ilen,lentmp
2055 character*3 out1file_text,ucase
2060 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2061 call getenv_loc("PREFIX",prefix)
2063 call getenv_loc("POT",pot)
2064 call getenv_loc("DIRTMP",tmpdir)
2065 call getenv_loc("CURDIR",curdir)
2066 call getenv_loc("OUT1FILE",out1file_text)
2067 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2068 out1file_text=ucase(out1file_text)
2069 if (out1file_text(1:1).eq."Y") then
2072 out1file=fg_rank.gt.0
2077 if (lentmp.gt.0) then
2078 write (*,'(80(1h!))')
2079 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2080 write (*,'(80(1h!))')
2081 write (*,*)"All output files will be on node /tmp directory."
2083 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2084 if (me.eq.king) then
2085 write (*,*) "The master node is ",nodename
2086 else if (fg_rank.eq.0) then
2087 write (*,*) "I am the CG slave node ",nodename
2089 write (*,*) "I am the FG slave node ",nodename
2092 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2093 lenpre = lentmp+lenpre+1
2095 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2096 C Get the names and open the input files
2097 #if defined(WINIFL) || defined(WINPGI)
2098 open(1,file=pref_orig(:ilen(pref_orig))//
2099 & '.inp',status='old',readonly,shared)
2100 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2101 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2102 C Get parameter filenames and open the parameter files.
2103 call getenv_loc('BONDPAR',bondname)
2104 open (ibond,file=bondname,status='old',readonly,shared)
2105 call getenv_loc('THETPAR',thetname)
2106 open (ithep,file=thetname,status='old',readonly,shared)
2107 call getenv_loc('ROTPAR',rotname)
2108 open (irotam,file=rotname,status='old',readonly,shared)
2109 call getenv_loc('TORPAR',torname)
2110 open (itorp,file=torname,status='old',readonly,shared)
2111 call getenv_loc('TORDPAR',tordname)
2112 open (itordp,file=tordname,status='old',readonly,shared)
2113 call getenv_loc('FOURIER',fouriername)
2114 open (ifourier,file=fouriername,status='old',readonly,shared)
2115 call getenv_loc('ELEPAR',elename)
2116 open (ielep,file=elename,status='old',readonly,shared)
2117 call getenv_loc('SIDEPAR',sidename)
2118 open (isidep,file=sidename,status='old',readonly,shared)
2119 call getenv_loc('LIPTRANPAR',liptranname)
2120 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2121 #elif (defined CRAY) || (defined AIX)
2122 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2124 c print *,"Processor",myrank," opened file 1"
2125 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2126 c print *,"Processor",myrank," opened file 9"
2127 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2128 C Get parameter filenames and open the parameter files.
2129 call getenv_loc('BONDPAR',bondname)
2130 open (ibond,file=bondname,status='old',action='read')
2131 c print *,"Processor",myrank," opened file IBOND"
2132 call getenv_loc('THETPAR',thetname)
2133 open (ithep,file=thetname,status='old',action='read')
2135 call getenv_loc('THETPARPDB',thetname_pdb)
2136 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2138 c print *,"Processor",myrank," opened file ITHEP"
2139 call getenv_loc('ROTPAR',rotname)
2140 open (irotam,file=rotname,status='old',action='read')
2142 call getenv_loc('ROTPARPDB',rotname_pdb)
2143 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2145 c print *,"Processor",myrank," opened file IROTAM"
2146 call getenv_loc('TORPAR',torname)
2147 open (itorp,file=torname,status='old',action='read')
2148 c print *,"Processor",myrank," opened file ITORP"
2149 call getenv_loc('TORDPAR',tordname)
2150 open (itordp,file=tordname,status='old',action='read')
2151 c print *,"Processor",myrank," opened file ITORDP"
2152 call getenv_loc('SCCORPAR',sccorname)
2153 open (isccor,file=sccorname,status='old',action='read')
2154 c print *,"Processor",myrank," opened file ISCCOR"
2155 call getenv_loc('FOURIER',fouriername)
2156 open (ifourier,file=fouriername,status='old',action='read')
2157 c print *,"Processor",myrank," opened file IFOURIER"
2158 call getenv_loc('ELEPAR',elename)
2159 open (ielep,file=elename,status='old',action='read')
2160 c print *,"Processor",myrank," opened file IELEP"
2161 call getenv_loc('SIDEPAR',sidename)
2162 open (isidep,file=sidename,status='old',action='read')
2163 call getenv_loc('LIPTRANPAR',liptranname)
2164 open (iliptranpar,file=liptranname,status='old',action='read')
2165 c print *,"Processor",myrank," opened file ISIDEP"
2166 c print *,"Processor",myrank," opened parameter files"
2168 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2169 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2170 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2171 C Get parameter filenames and open the parameter files.
2172 call getenv_loc('BONDPAR',bondname)
2173 open (ibond,file=bondname,status='old')
2174 call getenv_loc('THETPAR',thetname)
2175 open (ithep,file=thetname,status='old')
2177 call getenv_loc('THETPARPDB',thetname_pdb)
2178 open (ithep_pdb,file=thetname_pdb,status='old')
2180 call getenv_loc('ROTPAR',rotname)
2181 open (irotam,file=rotname,status='old')
2183 call getenv_loc('ROTPARPDB',rotname_pdb)
2184 open (irotam_pdb,file=rotname_pdb,status='old')
2186 call getenv_loc('TORPAR',torname)
2187 open (itorp,file=torname,status='old')
2188 call getenv_loc('TORDPAR',tordname)
2189 open (itordp,file=tordname,status='old')
2190 call getenv_loc('SCCORPAR',sccorname)
2191 open (isccor,file=sccorname,status='old')
2192 call getenv_loc('FOURIER',fouriername)
2193 open (ifourier,file=fouriername,status='old')
2194 call getenv_loc('ELEPAR',elename)
2195 open (ielep,file=elename,status='old')
2196 call getenv_loc('SIDEPAR',sidename)
2197 open (isidep,file=sidename,status='old')
2198 call getenv_loc('LIPTRANPAR',liptranname)
2199 open (iliptranpar,file=liptranname,status='old')
2201 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2203 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2204 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2205 C Get parameter filenames and open the parameter files.
2206 call getenv_loc('BONDPAR',bondname)
2207 open (ibond,file=bondname,status='old',readonly)
2208 call getenv_loc('THETPAR',thetname)
2209 open (ithep,file=thetname,status='old',readonly)
2210 call getenv_loc('ROTPAR',rotname)
2211 open (irotam,file=rotname,status='old',readonly)
2212 call getenv_loc('TORPAR',torname)
2213 open (itorp,file=torname,status='old',readonly)
2214 call getenv_loc('TORDPAR',tordname)
2215 open (itordp,file=tordname,status='old',readonly)
2216 call getenv_loc('SCCORPAR',sccorname)
2217 open (isccor,file=sccorname,status='old',readonly)
2219 call getenv_loc('THETPARPDB',thetname_pdb)
2220 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2222 call getenv_loc('FOURIER',fouriername)
2223 open (ifourier,file=fouriername,status='old',readonly)
2224 call getenv_loc('ELEPAR',elename)
2225 open (ielep,file=elename,status='old',readonly)
2226 call getenv_loc('SIDEPAR',sidename)
2227 open (isidep,file=sidename,status='old',readonly)
2228 call getenv_loc('LIPTRANPAR',liptranname)
2229 open (iliptranpar,file=liptranname,status='old',action='read')
2231 call getenv_loc('ROTPARPDB',rotname_pdb)
2232 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2236 call getenv_loc('TUBEPAR',tubename)
2237 #if defined(WINIFL) || defined(WINPGI)
2238 open (itube,file=tubename,status='old',readonly,shared)
2239 #elif (defined CRAY) || (defined AIX)
2240 open (itube,file=tubename,status='old',action='read')
2242 open (itube,file=tubename,status='old')
2244 open (itube,file=tubename,status='old',readonly)
2249 C 8/9/01 In the newest version SCp interaction constants are read from a file
2250 C Use -DOLDSCP to use hard-coded constants instead.
2252 call getenv_loc('SCPPAR',scpname)
2253 #if defined(WINIFL) || defined(WINPGI)
2254 open (iscpp,file=scpname,status='old',readonly,shared)
2255 #elif (defined CRAY) || (defined AIX)
2256 open (iscpp,file=scpname,status='old',action='read')
2258 open (iscpp,file=scpname,status='old')
2260 open (iscpp,file=scpname,status='old',readonly)
2263 call getenv_loc('PATTERN',patname)
2264 #if defined(WINIFL) || defined(WINPGI)
2265 open (icbase,file=patname,status='old',readonly,shared)
2266 #elif (defined CRAY) || (defined AIX)
2267 open (icbase,file=patname,status='old',action='read')
2269 open (icbase,file=patname,status='old')
2271 open (icbase,file=patname,status='old',readonly)
2274 C Open output file only for CG processes
2275 c print *,"Processor",myrank," fg_rank",fg_rank
2276 if (fg_rank.eq.0) then
2278 if (nodes.eq.1) then
2281 npos = dlog10(dfloat(nodes-1))+1
2283 if (npos.lt.3) npos=3
2284 write (liczba,'(i1)') npos
2285 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2287 write (liczba,form) me
2288 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2289 & liczba(:ilen(liczba))
2290 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2292 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2294 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2295 & liczba(:ilen(liczba))//'.mol2'
2296 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2297 & liczba(:ilen(liczba))//'.stat'
2299 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2300 & //liczba(:ilen(liczba))//'.stat')
2301 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2304 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2305 & liczba(:ilen(liczba))//'.const'
2310 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2311 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2312 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2313 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2314 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2316 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2318 rest2name=prefix(:ilen(prefix))//'.rst'
2320 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2323 #if defined(AIX) || defined(PGI) || defined(CRAY)
2324 if (me.eq.king .or. .not. out1file) then
2325 open(iout,file=outname,status='unknown')
2327 open(iout,file="/dev/null",status="unknown")
2331 if (fg_rank.gt.0) then
2332 write (liczba,'(i3.3)') myrank/nfgtasks
2333 write (ll,'(bz,i3.3)') fg_rank
2334 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2340 open(igeom,file=intname,status='unknown',position='append')
2341 open(ipdb,file=pdbname,status='unknown')
2342 open(imol2,file=mol2name,status='unknown')
2343 open(istat,file=statname,status='unknown',position='append')
2345 c1out open(iout,file=outname,status='unknown')
2348 if (me.eq.king .or. .not.out1file)
2349 & open(iout,file=outname,status='unknown')
2351 if (fg_rank.gt.0) then
2352 write (liczba,'(i3.3)') myrank/nfgtasks
2353 write (ll,'(bz,i3.3)') fg_rank
2354 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2359 open(igeom,file=intname,status='unknown',access='append')
2360 open(ipdb,file=pdbname,status='unknown')
2361 open(imol2,file=mol2name,status='unknown')
2362 open(istat,file=statname,status='unknown',access='append')
2364 c1out open(iout,file=outname,status='unknown')
2367 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2368 csa_seed=prefix(:lenpre)//'.CSA.seed'
2369 csa_history=prefix(:lenpre)//'.CSA.history'
2370 csa_bank=prefix(:lenpre)//'.CSA.bank'
2371 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2372 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2373 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2374 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2375 csa_int=prefix(:lenpre)//'.int'
2376 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2377 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2378 csa_in=prefix(:lenpre)//'.CSA.in'
2379 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2382 write (iout,'(80(1h-))')
2383 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2384 write (iout,'(80(1h-))')
2385 write (iout,*) "Input file : ",
2386 & pref_orig(:ilen(pref_orig))//'.inp'
2387 write (iout,*) "Output file : ",
2388 & outname(:ilen(outname))
2390 write (iout,*) "Sidechain potential file : ",
2391 & sidename(:ilen(sidename))
2393 write (iout,*) "SCp potential file : ",
2394 & scpname(:ilen(scpname))
2396 write (iout,*) "Electrostatic potential file : ",
2397 & elename(:ilen(elename))
2398 write (iout,*) "Cumulant coefficient file : ",
2399 & fouriername(:ilen(fouriername))
2400 write (iout,*) "Torsional parameter file : ",
2401 & torname(:ilen(torname))
2402 write (iout,*) "Double torsional parameter file : ",
2403 & tordname(:ilen(tordname))
2404 write (iout,*) "SCCOR parameter file : ",
2405 & sccorname(:ilen(sccorname))
2406 write (iout,*) "Bond & inertia constant file : ",
2407 & bondname(:ilen(bondname))
2408 write (iout,*) "Bending-torsion parameter file : ",
2409 & thetname(:ilen(thetname))
2410 write (iout,*) "Rotamer parameter file : ",
2411 & rotname(:ilen(rotname))
2412 write (iout,*) "Thetpdb parameter file : ",
2413 & thetname_pdb(:ilen(thetname_pdb))
2414 write (iout,*) "Threading database : ",
2415 & patname(:ilen(patname))
2417 &write (iout,*)" DIRTMP : ",
2419 write (iout,'(80(1h-))')
2423 c----------------------------------------------------------------------------
2424 subroutine card_concat(card)
2425 implicit real*8 (a-h,o-z)
2426 include 'DIMENSIONS'
2427 include 'COMMON.IOUNITS'
2429 character*80 karta,ucase
2431 read (inp,'(a)') karta
2434 do while (karta(80:80).eq.'&')
2435 card=card(:ilen(card)+1)//karta(:79)
2436 read (inp,'(a)') karta
2439 card=card(:ilen(card)+1)//karta
2442 c----------------------------------------------------------------------------------
2444 implicit real*8 (a-h,o-z)
2445 include 'DIMENSIONS'
2446 include 'COMMON.CHAIN'
2447 include 'COMMON.IOUNITS'
2449 open(irest2,file=rest2name,status='unknown')
2450 read(irest2,*) totT,EK,potE,totE,t_bath
2453 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2456 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2459 read (irest2,*) iset
2464 c---------------------------------------------------------------------------------
2465 subroutine read_fragments
2466 implicit real*8 (a-h,o-z)
2467 include 'DIMENSIONS'
2471 include 'COMMON.SETUP'
2472 include 'COMMON.CHAIN'
2473 include 'COMMON.IOUNITS'
2475 include 'COMMON.CONTROL'
2476 read(inp,*) nset,nfrag,npair,nfrag_back
2477 loc_qlike=(nfrag_back.lt.0)
2478 nfrag_back=iabs(nfrag_back)
2479 if(me.eq.king.or..not.out1file)
2480 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2481 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2483 read(inp,*) mset(iset)
2485 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2487 if(me.eq.king.or..not.out1file)
2488 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2489 & ifrag(2,i,iset), qinfrag(i,iset)
2492 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2494 if(me.eq.king.or..not.out1file)
2495 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2496 & ipair(2,i,iset), qinpair(i,iset)
2500 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2501 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2502 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2503 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2504 if(me.eq.king.or..not.out1file)
2505 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2506 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2507 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2508 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2512 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2513 & wfrag_back(3,i,iset),
2514 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2515 if(me.eq.king.or..not.out1file)
2516 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2517 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2523 C---------------------------------------------------------------------------
2524 subroutine read_afminp
2525 implicit real*8 (a-h,o-z)
2526 include 'DIMENSIONS'
2530 include 'COMMON.SETUP'
2531 include 'COMMON.CONTROL'
2532 include 'COMMON.CHAIN'
2533 include 'COMMON.IOUNITS'
2534 include 'COMMON.SBRIDGE'
2535 character*320 afmcard
2537 call card_concat(afmcard)
2538 call readi(afmcard,"BEG",afmbeg,0)
2539 call readi(afmcard,"END",afmend,0)
2540 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2541 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2542 print *,'FORCE=' ,forceAFMconst
2543 CCCC NOW PROPERTIES FOR AFM
2546 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2548 distafminit=dsqrt(distafminit)
2549 print *,'initdist',distafminit
2552 c-------------------------------------------------------------------------------
2553 subroutine read_saxs_constr
2554 implicit real*8 (a-h,o-z)
2555 include 'DIMENSIONS'
2559 include 'COMMON.SETUP'
2560 include 'COMMON.CONTROL'
2561 include 'COMMON.CHAIN'
2562 include 'COMMON.IOUNITS'
2563 include 'COMMON.SBRIDGE'
2564 double precision cm(3)
2566 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2568 if (saxs_mode.eq.0) then
2569 c SAXS distance distribution
2571 read(inp,*) distsaxs(i),Psaxs(i)
2575 Cnorm = Cnorm + Psaxs(i)
2577 write (iout,*) "Cnorm",Cnorm
2579 Psaxs(i)=Psaxs(i)/Cnorm
2581 write (iout,*) "Normalized distance distribution from SAXS"
2583 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2587 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2589 write (iout,*) "Wsaxs0",Wsaxs0
2593 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2600 cm(j)=cm(j)+Csaxs(j,i)
2608 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2611 write (iout,*) "SAXS sphere coordinates"
2613 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2619 c-------------------------------------------------------------------------------
2620 subroutine read_dist_constr
2621 implicit real*8 (a-h,o-z)
2622 include 'DIMENSIONS'
2626 include 'COMMON.SETUP'
2627 include 'COMMON.CONTROL'
2628 include 'COMMON.CHAIN'
2629 include 'COMMON.IOUNITS'
2630 include 'COMMON.SBRIDGE'
2631 include 'COMMON.INTERACT'
2632 integer ifrag_(2,100),ipair_(2,1000)
2633 double precision wfrag_(100),wpair_(1000)
2634 character*5000 controlcard
2635 logical normalize,next
2637 double precision scal_bfac
2638 double precision xlink(4,0:4) /
2640 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2641 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2642 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2643 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2644 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2645 c print *, "WCHODZE"
2646 write (iout,*) "Calling read_dist_constr"
2647 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2649 restr_on_coord=.false.
2658 call card_concat(controlcard)
2659 next = index(controlcard,"NEXT").gt.0
2660 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2661 write (iout,*) "restr_type",restr_type
2662 call readi(controlcard,"NFRAG",nfrag_,0)
2663 call readi(controlcard,"NFRAG",nfrag_,0)
2664 call readi(controlcard,"NPAIR",npair_,0)
2665 call readi(controlcard,"NDIST",ndist_,0)
2666 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2667 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2668 if (restr_type.eq.10)
2669 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2670 if (restr_type.eq.12)
2671 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2672 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2673 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2674 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2675 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2676 normalize = index(controlcard,"NORMALIZE").gt.0
2677 write (iout,*) "WBOLTZD",wboltzd
2678 write (iout,*) "SCAL_PEAK",scal_peak
2679 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2680 write (iout,*) "IFRAG"
2682 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2684 write (iout,*) "IPAIR"
2686 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2688 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2690 & "Distance restraints as generated from reference structure"
2692 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2693 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2694 & ifrag_(2,i)=nstart_sup+nsup-1
2695 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2697 if (wfrag_(i).eq.0.0d0) cycle
2698 do j=ifrag_(1,i),ifrag_(2,i)-1
2699 do k=j+1,ifrag_(2,i)
2700 c write (iout,*) "j",j," k",k
2702 if (restr_type.eq.1) then
2708 forcon(nhpb)=wfrag_(i)
2709 else if (constr_dist.eq.2) then
2710 if (ddjk.le.dist_cut) then
2716 forcon(nhpb)=wfrag_(i)
2718 else if (restr_type.eq.3) then
2724 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2727 if (.not.out1file .or. me.eq.king)
2728 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2729 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2731 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2732 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2738 if (wpair_(i).eq.0.0d0) cycle
2746 do j=ifrag_(1,ii),ifrag_(2,ii)
2747 do k=ifrag_(1,jj),ifrag_(2,jj)
2749 if (restr_type.eq.1) then
2755 forcon(nhpb)=wpair_(i)
2756 else if (constr_dist.eq.2) then
2757 if (ddjk.le.dist_cut) then
2763 forcon(nhpb)=wpair_(i)
2765 else if (restr_type.eq.3) then
2771 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2774 if (.not.out1file .or. me.eq.king)
2775 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2776 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2778 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2779 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2786 write (iout,*) "Distance restraints as read from input"
2788 if (restr_type.eq.12) then
2789 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2790 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2791 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2792 & fordepth_peak(nhpb_peak+1),npeak
2793 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2794 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2795 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2796 c & fordepth_peak(nhpb_peak+1),npeak
2797 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2798 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2799 nhpb_peak=nhpb_peak+1
2800 irestr_type_peak(nhpb_peak)=12
2801 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2804 if (.not.out1file .or. me.eq.king)
2805 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2806 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2807 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2808 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2809 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2811 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2812 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2813 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2814 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2815 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2817 if (ibecarb_peak(nhpb_peak).eq.3) then
2818 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2819 else if (ibecarb_peak(nhpb_peak).eq.2) then
2820 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2821 else if (ibecarb_peak(nhpb_peak).eq.1) then
2822 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2823 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2825 else if (restr_type.eq.11) then
2826 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2827 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2828 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2829 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2831 irestr_type(nhpb)=11
2833 if (.not.out1file .or. me.eq.king)
2834 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2835 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2836 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2838 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2839 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2840 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2842 c if (ibecarb(nhpb).gt.0) then
2843 c ihpb(nhpb)=ihpb(nhpb)+nres
2844 c jhpb(nhpb)=jhpb(nhpb)+nres
2846 if (ibecarb(nhpb).eq.3) then
2847 ihpb(nhpb)=ihpb(nhpb)+nres
2848 else if (ibecarb(nhpb).eq.2) then
2849 ihpb(nhpb)=ihpb(nhpb)+nres
2850 else if (ibecarb(nhpb).eq.1) then
2851 ihpb(nhpb)=ihpb(nhpb)+nres
2852 jhpb(nhpb)=jhpb(nhpb)+nres
2854 else if (restr_type.eq.10) then
2855 c Cross-lonk Markov-like potential
2856 call card_concat(controlcard)
2857 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2858 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2860 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2861 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2862 if (index(controlcard,"ZL").gt.0) then
2864 else if (index(controlcard,"ADH").gt.0) then
2866 else if (index(controlcard,"PDH").gt.0) then
2868 else if (index(controlcard,"DSS").gt.0) then
2873 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2874 & xlink(1,link_type))
2875 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2876 & xlink(2,link_type))
2877 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2878 & xlink(3,link_type))
2879 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2880 & xlink(4,link_type))
2881 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2882 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2883 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2884 if (forcon(nhpb+1).le.0.0d0 .or.
2885 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2887 irestr_type(nhpb)=10
2888 if (ibecarb(nhpb).eq.3) then
2889 jhpb(nhpb)=jhpb(nhpb)+nres
2890 else if (ibecarb(nhpb).eq.2) then
2891 ihpb(nhpb)=ihpb(nhpb)+nres
2892 else if (ibecarb(nhpb).eq.1) then
2893 ihpb(nhpb)=ihpb(nhpb)+nres
2894 jhpb(nhpb)=jhpb(nhpb)+nres
2897 if (.not.out1file .or. me.eq.king)
2898 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2899 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2900 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2903 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2904 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2905 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2910 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2911 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2912 if (forcon(nhpb+1).gt.0.0d0) then
2914 if (dhpb1(nhpb).eq.0.0d0) then
2919 if (ibecarb(nhpb).eq.3) then
2920 jhpb(nhpb)=jhpb(nhpb)+nres
2921 else if (ibecarb(nhpb).eq.2) then
2922 ihpb(nhpb)=ihpb(nhpb)+nres
2923 else if (ibecarb(nhpb).eq.1) then
2924 ihpb(nhpb)=ihpb(nhpb)+nres
2925 jhpb(nhpb)=jhpb(nhpb)+nres
2927 if (dhpb(nhpb).eq.0.0d0)
2928 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2931 if (.not.out1file .or. me.eq.king)
2932 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2933 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2935 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2936 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2939 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2940 C if (forcon(nhpb+1).gt.0.0d0) then
2942 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2945 if (restr_type.eq.4) then
2946 write (iout,*) "The BFAC array"
2948 write (iout,'(i5,f10.5)') i,bfac(i)
2951 if (itype(i).eq.ntyp1) cycle
2953 if (itype(j).eq.ntyp1) cycle
2954 if (itype(i).eq.10) then
2959 if (itype(j).eq.10) then
2969 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2972 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2973 ihpb(nhpb)=i+nres*ii
2974 jhpb(nhpb)=j+nres*jj
2975 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2977 if (.not.out1file .or. me.eq.king) then
2978 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2979 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2980 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2984 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2985 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2986 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2996 if (restr_type.eq.5) then
2997 restr_on_coord=.true.
2999 if (itype(i).eq.ntyp1) cycle
3000 bfac(i)=(scal_bfac/bfac(i))**2
3009 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3010 & fordepthmax=fordepth(i)
3013 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3018 c-------------------------------------------------------------------------------
3019 subroutine read_constr_homology
3021 include 'DIMENSIONS'
3025 include 'COMMON.SETUP'
3026 include 'COMMON.CONTROL'
3027 include 'COMMON.CHAIN'
3028 include 'COMMON.IOUNITS'
3030 include 'COMMON.GEO'
3031 include 'COMMON.INTERACT'
3032 include 'COMMON.NAMES'
3034 c For new homol impl
3036 include 'COMMON.VAR'
3039 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3041 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3042 c & sigma_odl_temp(maxres,maxres,max_template)
3044 character*24 model_ki_dist, model_ki_angle
3045 character*500 controlcard
3046 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3051 c FP - Nov. 2014 Temporary specifications for new vars
3053 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3055 double precision, dimension (max_template,maxres) :: rescore
3056 double precision, dimension (max_template,maxres) :: rescore2
3057 double precision, dimension (max_template,maxres) :: rescore3
3058 character*24 pdbfile,tpl_k_rescore
3059 c -----------------------------------------------------------------
3060 c Reading multiple PDB ref structures and calculation of retraints
3061 c not using pre-computed ones stored in files model_ki_{dist,angle}
3063 c -----------------------------------------------------------------
3066 c Alternative: reading from input
3067 call card_concat(controlcard)
3068 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3069 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3070 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3071 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3072 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3073 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3074 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3075 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3076 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3077 if(.not.read2sigma.and.start_from_model) then
3078 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3079 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3080 start_from_model=.false.
3082 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3083 & write(iout,*) 'START_FROM_MODELS is ON'
3084 if(start_from_model .and. rest) then
3085 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3086 write(iout,*) 'START_FROM_MODELS is OFF'
3087 write(iout,*) 'remove restart keyword from input'
3090 if (homol_nset.gt.1)then
3091 call card_concat(controlcard)
3092 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3093 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3094 write(iout,*) "iset homology_weight "
3096 write(iout,*) i,waga_homology(i)
3099 iset=mod(kolor,homol_nset)+1
3102 waga_homology(1)=1.0
3105 cd write (iout,*) "nnt",nnt," nct",nct
3112 c write(iout,*) 'nnt=',nnt,'nct=',nct
3115 do k=1,constr_homology
3128 if (read_homol_frag) then
3129 call read_klapaucjusz
3132 do k=1,constr_homology
3134 read(inp,'(a)') pdbfile
3135 if(me.eq.king .or. .not. out1file)
3136 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3137 & pdbfile(:ilen(pdbfile))
3138 open(ipdbin,file=pdbfile,status='old',err=33)
3140 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3141 & pdbfile(:ilen(pdbfile))
3144 c print *,'Begin reading pdb data'
3146 c Files containing res sim or local scores (former containing sigmas)
3149 write(kic2,'(bz,i2.2)') k
3151 tpl_k_rescore="template"//kic2//".sco"
3154 if (read2sigma) then
3155 call readpdb_template(k)
3160 c Distance restraints
3163 C Copy the coordinates from reference coordinates (?)
3167 c write (iout,*) "c(",j,i,") =",c(j,i)
3171 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3173 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3174 open (ientin,file=tpl_k_rescore,status='old')
3175 if (nnt.gt.1) rescore(k,1)=0.0d0
3176 do irec=nnt,nct ! loop for reading res sim
3177 if (read2sigma) then
3178 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3179 & rescore3_tmp,idomain_tmp
3181 idomain(k,i_tmp)=idomain_tmp
3182 rescore(k,i_tmp)=rescore_tmp
3183 rescore2(k,i_tmp)=rescore2_tmp
3184 rescore3(k,i_tmp)=rescore3_tmp
3185 if (.not. out1file .or. me.eq.king)
3186 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3187 & i_tmp,rescore2_tmp,rescore_tmp,
3188 & rescore3_tmp,idomain_tmp
3191 read (ientin,*,end=1401) rescore_tmp
3193 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3194 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3195 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3200 if (waga_dist.ne.0.0d0) then
3208 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3209 c write (iout,*) k,i,j,distal,dist2_cut
3211 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3212 & .and. distal.le.dist2_cut ) then
3218 c write (iout,*) "k",k
3219 c write (iout,*) "i",i," j",j," constr_homology",
3224 if (read2sigma) then
3227 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3229 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3230 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3231 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3233 if (odl(k,ii).le.dist_cut) then
3234 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3237 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3238 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3240 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3241 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3245 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3248 l_homo(k,ii)=.false.
3255 c Theta, dihedral and SC retraints
3257 if (waga_angle.gt.0.0d0) then
3258 c open (ientin,file=tpl_k_sigma_dih,status='old')
3259 c do irec=1,maxres-3 ! loop for reading sigma_dih
3260 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3261 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3262 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3263 c & sigma_dih(k,i+nnt-1)
3268 if (idomain(k,i).eq.0) then
3272 dih(k,i)=phiref(i) ! right?
3273 c read (ientin,*) sigma_dih(k,i) ! original variant
3274 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3275 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3276 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3277 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3278 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3280 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3281 & rescore(k,i-2)+rescore(k,i-3))/4.0
3282 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3283 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3284 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3285 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3286 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3287 c Instead of res sim other local measure of b/b str reliability possible
3288 if (sigma_dih(k,i).ne.0)
3289 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3290 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3295 if (waga_theta.gt.0.0d0) then
3296 c open (ientin,file=tpl_k_sigma_theta,status='old')
3297 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3298 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3299 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3300 c & sigma_theta(k,i+nnt-1)
3305 do i = nnt+2,nct ! right? without parallel.
3306 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3307 c do i=ithet_start,ithet_end ! with FG parallel.
3308 if (idomain(k,i).eq.0) then
3309 sigma_theta(k,i)=0.0
3312 thetatpl(k,i)=thetaref(i)
3313 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3314 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3315 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3316 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3317 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3318 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3319 & rescore(k,i-2))/3.0
3320 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3321 if (sigma_theta(k,i).ne.0)
3322 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3324 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3325 c rescore(k,i-2) ! right expression ?
3326 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3330 if (waga_d.gt.0.0d0) then
3331 c open (ientin,file=tpl_k_sigma_d,status='old')
3332 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3333 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3334 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3335 c & sigma_d(k,i+nnt-1)
3339 do i = nnt,nct ! right? without parallel.
3340 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3341 c do i=loc_start,loc_end ! with FG parallel.
3342 if (itype(i).eq.10) cycle
3343 if (idomain(k,i).eq.0 ) then
3350 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3351 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3352 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3353 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3354 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3355 if (sigma_d(k,i).ne.0)
3356 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3358 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3359 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3360 c read (ientin,*) sigma_d(k,i) ! 1st variant
3365 c remove distance restraints not used in any model from the list
3366 c shift data in all arrays
3368 if (waga_dist.ne.0.0d0) then
3374 if (ii_in_use(ii).eq.0.and.liiflag) then
3378 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3379 & .not.liiflag.and.ii.eq.lim_odl) then
3380 if (ii.eq.lim_odl) then
3381 iishift=ii-iistart+1
3386 do ki=iistart,lim_odl-iishift
3387 ires_homo(ki)=ires_homo(ki+iishift)
3388 jres_homo(ki)=jres_homo(ki+iishift)
3389 ii_in_use(ki)=ii_in_use(ki+iishift)
3390 do k=1,constr_homology
3391 odl(k,ki)=odl(k,ki+iishift)
3392 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3393 l_homo(k,ki)=l_homo(k,ki+iishift)
3397 lim_odl=lim_odl-iishift
3403 endif ! .not. klapaucjusz
3405 if (constr_homology.gt.0) call homology_partition
3406 if (constr_homology.gt.0) call init_int_table
3407 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3408 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3412 if (.not.out_template_restr) return
3413 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3414 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3415 write (iout,*) "Distance restraints from templates"
3417 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3418 & ii,ires_homo(ii),jres_homo(ii),
3419 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3420 & ki=1,constr_homology)
3422 write (iout,*) "Dihedral angle restraints from templates"
3424 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3425 & (rad2deg*dih(ki,i),
3426 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3428 write (iout,*) "Virtual-bond angle restraints from templates"
3430 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3431 & (rad2deg*thetatpl(ki,i),
3432 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3434 write (iout,*) "SC restraints from templates"
3436 write(iout,'(i5,100(4f8.2,4x))') i,
3437 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3438 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3441 c -----------------------------------------------------------------
3444 c----------------------------------------------------------------------
3446 subroutine flush(iu)
3451 subroutine flush(iu)
3456 c------------------------------------------------------------------------------
3457 subroutine copy_to_tmp(source)
3458 include "DIMENSIONS"
3459 include "COMMON.IOUNITS"
3460 character*(*) source
3461 character* 256 tmpfile
3465 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3466 inquire(file=tmpfile,exist=ex)
3468 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3469 & " to temporary directory..."
3470 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3471 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3475 c------------------------------------------------------------------------------
3476 subroutine move_from_tmp(source)
3477 include "DIMENSIONS"
3478 include "COMMON.IOUNITS"
3479 character*(*) source
3482 write (*,*) "Moving ",source(:ilen(source)),
3483 & " from temporary directory to working directory"
3484 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3485 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3488 c------------------------------------------------------------------------------
3489 subroutine random_init(seed)
3491 C Initialize random number generator
3493 implicit real*8 (a-h,o-z)
3494 include 'DIMENSIONS'
3497 logical OKRandom, prng_restart
3499 integer iseed_array(4)
3501 include 'COMMON.IOUNITS'
3502 include 'COMMON.TIME1'
3503 include 'COMMON.THREAD'
3504 include 'COMMON.SBRIDGE'
3505 include 'COMMON.CONTROL'
3506 include 'COMMON.MCM'
3507 include 'COMMON.MAP'
3508 include 'COMMON.HEADER'
3509 include 'COMMON.CSA'
3510 include 'COMMON.CHAIN'
3511 include 'COMMON.MUCA'
3513 include 'COMMON.FFIELD'
3514 include 'COMMON.SETUP'
3515 iseed=-dint(dabs(seed))
3516 if (iseed.eq.0) then
3517 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3518 & 'Random seed undefined. The program will stop.'
3519 write (*,'(/80(1h*)/20x,a/80(1h*))')
3520 & 'Random seed undefined. The program will stop.'
3522 call mpi_finalize(mpi_comm_world,ierr)
3524 stop 'Bad random seed.'
3527 if (fg_rank.eq.0) then
3531 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3532 OKRandom = prng_restart(me,iseed)
3535 tmp=65536.0d0**(4-i)
3536 iseed_array(i) = dint(seed/tmp)
3537 seed=seed-iseed_array(i)*tmp
3540 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3541 & (iseed_array(i),i=1,4)
3542 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3543 & (iseed_array(i),i=1,4)
3544 OKRandom = prng_restart(me,iseed_array)
3547 c r1 = prng_next(me)
3548 r1=ran_number(0.0D0,1.0D0)
3550 & write (iout,*) 'ran_num',r1
3551 if (r1.lt.0.0d0) OKRandom=.false.
3553 if (.not.OKRandom) then
3554 write (iout,*) 'PRNG IS NOT WORKING!!!'
3555 print *,'PRNG IS NOT WORKING!!!'
3558 call mpi_abort(mpi_comm_world,error_msg,ierr)
3561 write (iout,*) 'too many processors for parallel prng'
3562 write (*,*) 'too many processors for parallel prng'
3570 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3574 c----------------------------------------------------------------------
3575 subroutine read_klapaucjusz
3577 include 'DIMENSIONS'
3581 include 'COMMON.SETUP'
3582 include 'COMMON.CONTROL'
3583 include 'COMMON.CHAIN'
3584 include 'COMMON.IOUNITS'
3586 include 'COMMON.GEO'
3587 include 'COMMON.INTERACT'
3588 include 'COMMON.NAMES'
3589 character*256 fragfile
3590 integer ninclust(maxclust),inclust(max_template,maxclust),
3591 & nresclust(maxclust),iresclust(maxres,maxclust)
3594 character*24 model_ki_dist, model_ki_angle
3595 character*500 controlcard
3596 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3597 logical lprn /.true./
3603 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3604 double precision, dimension (max_template,maxres) :: rescore
3605 double precision, dimension (max_template,maxres) :: rescore2
3606 character*24 pdbfile,tpl_k_rescore
3609 c For new homol impl
3611 include 'COMMON.VAR'
3613 call getenv("FRAGFILE",fragfile)
3614 open(ientin,file=fragfile,status="old",err=10)
3615 read(ientin,*) constr_homology,nclust
3621 do k=1,constr_homology
3622 read(ientin,'(a)') pdbfile
3623 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3624 & pdbfile(:ilen(pdbfile))
3625 open(ipdbin,file=pdbfile,status='old',err=33)
3627 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3628 & pdbfile(:ilen(pdbfile))
3632 call readpdb_template(k)
3640 read(ientin,*) ninclust(i),nresclust(i)
3641 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3642 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3645 c Loop over clusters
3648 do ll = 1,ninclust(l)
3656 idomain(k,iresclust(i,l)+1) = 1
3658 idomain(k,iresclust(i,l)) = 1
3662 c Distance restraints
3665 C Copy the coordinates from reference coordinates (?)
3669 c write (iout,*) "c(",j,i,") =",c(j,i)
3672 call int_from_cart(.true.,.false.)
3673 call sc_loc_geom(.true.)
3675 thetaref(i)=theta(i)
3678 if (waga_dist.ne.0.0d0) then
3686 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3687 c write (iout,*) k,i,j,distal,dist2_cut
3689 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3690 & .and. distal.le.dist2_cut ) then
3696 c write (iout,*) "k",k
3697 c write (iout,*) "i",i," j",j," constr_homology",
3702 if (read2sigma) then
3705 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3707 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3708 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3709 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3711 if (odl(k,ii).le.dist_cut) then
3712 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3715 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3716 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3718 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3719 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3723 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3726 c l_homo(k,ii)=.false.
3733 c Theta, dihedral and SC retraints
3735 if (waga_angle.gt.0.0d0) then
3737 if (idomain(k,i).eq.0) then
3738 c sigma_dih(k,i)=0.0
3742 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3743 & rescore(k,i-2)+rescore(k,i-3))/4.0
3744 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3745 c & " sigma_dihed",sigma_dih(k,i)
3746 if (sigma_dih(k,i).ne.0)
3747 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3752 if (waga_theta.gt.0.0d0) then
3754 if (idomain(k,i).eq.0) then
3755 c sigma_theta(k,i)=0.0
3758 thetatpl(k,i)=thetaref(i)
3759 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3760 & rescore(k,i-2))/3.0
3761 if (sigma_theta(k,i).ne.0)
3762 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3766 if (waga_d.gt.0.0d0) then
3768 if (itype(i).eq.10) cycle
3769 if (idomain(k,i).eq.0 ) then
3776 sigma_d(k,i)=rescore(k,i)
3777 if (sigma_d(k,i).ne.0)
3778 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3779 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3785 c remove distance restraints not used in any model from the list
3786 c shift data in all arrays
3788 if (waga_dist.ne.0.0d0) then
3794 if (ii_in_use(ii).eq.0.and.liiflag) then
3798 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3799 & .not.liiflag.and.ii.eq.lim_odl) then
3800 if (ii.eq.lim_odl) then
3801 iishift=ii-iistart+1
3806 do ki=iistart,lim_odl-iishift
3807 ires_homo(ki)=ires_homo(ki+iishift)
3808 jres_homo(ki)=jres_homo(ki+iishift)
3809 ii_in_use(ki)=ii_in_use(ki+iishift)
3810 do k=1,constr_homology
3811 odl(k,ki)=odl(k,ki+iishift)
3812 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3813 l_homo(k,ki)=l_homo(k,ki+iishift)
3817 lim_odl=lim_odl-iishift
3824 10 stop "Error in fragment file"