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)
724 C Read PDB structure if applicable
726 if (indpdb.gt.0 .or. pdbref) then
727 read(inp,'(a)') pdbfile
728 if(me.eq.king.or..not.out1file)
729 & write (iout,'(2a)') 'PDB data will be read from file ',
730 & pdbfile(:ilen(pdbfile))
731 open(ipdbin,file=pdbfile,status='old',err=33)
733 33 write (iout,'(a)') 'Error opening PDB file.'
744 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
745 & (crefjlee(j,i+nres),j=1,3)
748 if(me.eq.king.or..not.out1file)
749 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
750 & ' nstart_sup=',nstart_sup
752 itype_pdb(i)=itype(i)
756 nct=nstart_sup+nsup-1
757 call contact(.false.,ncont_ref,icont_ref,co)
760 if(me.eq.king.or..not.out1file)
761 & write(iout,*)'Adding sidechains'
765 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
768 do while (fail.and.nsi.le.maxsi)
769 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
772 if(fail) write(iout,*)'Adding sidechain failed for res ',
773 & i,' after ',nsi,' trials'
778 if (indpdb.eq.0) then
779 C Read sequence if not taken from the pdb file.
781 c print *,'nres=',nres
782 if (iscode.gt.0) then
783 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
785 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
787 C Convert sequence to numeric code
789 itype(i)=rescode(i,sequence(i),iscode)
791 C Assign initial virtual bond lengths
797 vbld(i+nres)=dsc(iabs(itype(i)))
798 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
799 c write (iout,*) "i",i," itype",itype(i),
800 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
804 c print '(20i4)',(itype(i),i=1,nres)
807 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
809 if (itype(i).eq.ntyp1) then
813 else if (iabs(itype(i+1)).ne.20) then
815 else if (iabs(itype(i)).ne.20) then
822 if(me.eq.king.or..not.out1file)then
823 write (iout,*) "ITEL"
825 write (iout,*) i,itype(i),itel(i)
827 print *,'Call Read_Bridge.'
831 cd print *,'NNT=',NNT,' NCT=',NCT
832 call seq2chains(nres,itype,nchain,chain_length,chain_border,
834 write(iout,*) "nres",nres," nchain",nchain
836 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
839 call chain_symmetry(nchain,nres,itype,chain_border,
840 & chain_length,npermchain,tabpermchain)
842 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
845 write(iout,*) "residue permutations"
847 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
849 if (itype(1).eq.ntyp1) nnt=2
850 if (itype(nres).eq.ntyp1) nct=nct-1
852 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
853 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
855 print*, 'init_dfa_vars finished!'
857 print*, 'read_dfa_info finished!'
861 if(me.eq.king.or..not.out1file)
862 & write (iout,'(a,i3)') 'nsup=',nsup
864 if (nsup.le.(nct-nnt+1)) then
865 do i=0,nct-nnt+1-nsup
866 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
872 & 'Error - sequences to be superposed do not match.'
875 do i=0,nsup-(nct-nnt+1)
876 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
878 nstart_sup=nstart_sup+i
884 & 'Error - sequences to be superposed do not match.'
887 if (nsup.eq.0) nsup=nct-nnt
888 if (nstart_sup.eq.0) nstart_sup=nnt
889 if (nstart_seq.eq.0) nstart_seq=nnt
890 if(me.eq.king.or..not.out1file)
891 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
892 & ' nstart_seq=',nstart_seq
895 C 8/13/98 Set limits to generating the dihedral angles
900 read (inp,*) ndih_constr
901 write (iout,*) "ndih_constr",ndih_constr
902 if (ndih_constr.gt.0) then
905 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
908 if(me.eq.king.or..not.out1file)then
911 & 'There are',ndih_constr,' restraints on gamma angles.'
913 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
920 phi0(i)=deg2rad*phi0(i)
921 drange(i)=deg2rad*drange(i)
923 C if(me.eq.king.or..not.out1file)
924 C & write (iout,*) 'FTORS',ftors
927 phibound(1,ii) = phi0(i)-drange(i)
928 phibound(2,ii) = phi0(i)+drange(i)
931 if (me.eq.king .or. .not.out1file) then
933 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
935 write (iout,'(a3,i5,2f10.1)')
936 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
941 else if (ndih_constr.lt.0) then
943 call card_concat(weightcard)
944 call reada(weightcard,"PHIHEL",phihel,50.0D0)
945 call reada(weightcard,"PHIBET",phibet,180.0D0)
946 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
947 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
948 call reada(weightcard,"WDIHC",wdihc,0.591D0)
949 write (iout,*) "Weight of dihedral angle restraints",wdihc
950 read(inp,'(9x,3f7.3)')
951 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
952 write (iout,*) "The secprob array"
954 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
958 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
959 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
960 ndih_constr=ndih_constr+1
961 idih_constr(ndih_constr)=i
964 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
965 sumv=sumv+vpsipred(j,ndih_constr)
968 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
970 phibound(1,ndih_constr)=phihel*deg2rad
971 phibound(2,ndih_constr)=phibet*deg2rad
972 sdihed(1,ndih_constr)=sigmahel*deg2rad
973 sdihed(2,ndih_constr)=sigmabet*deg2rad
977 if(me.eq.king.or..not.out1file)then
980 & 'There are',ndih_constr,
981 & ' bimodal restraints on gamma angles.'
983 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
984 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
985 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
986 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
987 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
988 & (vpsipred(j,i),j=1,3)
995 C first setting the theta boundaries to 0 to pi
996 C this mean that there is no energy penalty for any angle occuring this can be applied
997 C for generate random conformation but is not implemented in this way
1000 C thetabound(2,i)=pi
1002 C begin reading theta constrains this is quartic constrains allowing to
1003 C have smooth second derivative
1004 if (with_theta_constr) then
1005 C with_theta_constr is keyword allowing for occurance of theta constrains
1006 read (inp,*) ntheta_constr
1007 C ntheta_constr is the number of theta constrains
1008 if (ntheta_constr.gt.0) then
1009 C read (inp,*) ftors
1010 read (inp,*) (itheta_constr(i),theta_constr0(i),
1011 & theta_drange(i),for_thet_constr(i),
1012 & i=1,ntheta_constr)
1013 C the above code reads from 1 to ntheta_constr
1014 C itheta_constr(i) residue i for which is theta_constr
1015 C theta_constr0 the global minimum value
1016 C theta_drange is range for which there is no energy penalty
1017 C for_thet_constr is the force constant for quartic energy penalty
1019 if(me.eq.king.or..not.out1file)then
1021 & 'There are',ntheta_constr,' constraints on phi angles.'
1022 do i=1,ntheta_constr
1023 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1025 & for_thet_constr(i)
1028 do i=1,ntheta_constr
1029 theta_constr0(i)=deg2rad*theta_constr0(i)
1030 theta_drange(i)=deg2rad*theta_drange(i)
1032 C if(me.eq.king.or..not.out1file)
1033 C & write (iout,*) 'FTORS',ftors
1034 C do i=1,ntheta_constr
1035 C ii = itheta_constr(i)
1036 C thetabound(1,ii) = phi0(i)-drange(i)
1037 C thetabound(2,ii) = phi0(i)+drange(i)
1039 endif ! ntheta_constr.gt.0
1040 endif! with_theta_constr
1042 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1043 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1044 c--- Zscore rms -------
1045 if (nz_start.eq.0) nz_start=nnt
1046 if (nz_end.eq.0 .and. nsup.gt.0) then
1048 else if (nz_end.eq.0) then
1051 if(me.eq.king.or..not.out1file)then
1052 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1053 write (iout,*) 'IZ_SC=',iz_sc
1055 c----------------------
1058 if (.not.pdbref) then
1059 call read_angles(inp,*38)
1061 38 write (iout,'(a)') 'Error reading reference structure.'
1063 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1064 stop 'Error reading reference structure'
1066 39 call chainbuild_extconf
1068 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1077 call contact(.true.,ncont_ref,icont_ref,co)
1081 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1083 if (constr_dist.gt.0) call read_dist_constr
1084 write (iout,*) "After read_dist_constr nhpb",nhpb
1085 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1087 call NMRpeak_partition
1088 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1090 if(me.eq.king.or..not.out1file)
1091 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1094 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1096 if(me.eq.king.or..not.out1file)
1097 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1098 & icont_ref(1,i),' ',
1099 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1102 write (iout,*) "calling read_saxs_consrtr",nsaxs
1103 if (nsaxs.gt.0) call read_saxs_constr
1105 if (constr_homology.gt.0) then
1106 call read_constr_homology
1107 if (indpdb.gt.0 .or. pdbref) then
1110 c(j,i)=crefjlee(j,i)
1111 cref(j,i)=crefjlee(j,i)
1116 write (iout,*) "sc_loc_geom: Array C"
1118 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1119 & (c(j,i+nres),j=1,3)
1121 write (iout,*) "Array Cref"
1123 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1124 & (cref(j,i+nres),j=1,3)
1127 call int_from_cart1(.false.)
1128 call sc_loc_geom(.false.)
1130 thetaref(i)=theta(i)
1135 dc(j,i)=c(j,i+1)-c(j,i)
1136 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1141 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1142 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1151 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1152 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1153 & modecalc.ne.10) then
1154 C If input structure hasn't been supplied from the PDB file read or generate
1156 if (iranconf.eq.0 .and. .not. extconf) then
1157 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1158 & write (iout,'(a)') 'Initial geometry will be read in.'
1160 read(inp,'(8f10.5)',end=36,err=36)
1161 & ((c(l,k),l=1,3),k=1,nres),
1162 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1163 write (iout,*) "Exit READ_CART"
1164 c write (iout,'(8f10.5)')
1165 c & ((c(l,k),l=1,3),k=1,nres),
1166 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1170 dc(j,i)=c(j,i+1)-c(j,i)
1171 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1175 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1177 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1178 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1183 c dc_norm(j,i+nres)=0.0d0
1187 call int_from_cart1(.true.)
1188 write (iout,*) "Finish INT_TO_CART"
1189 c write (iout,*) "DC"
1191 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1192 c & (dc(j,i+nres),j=1,3)
1198 call read_angles(inp,*36)
1199 call chainbuild_extconf
1202 36 write (iout,'(a)') 'Error reading angle file.'
1204 call mpi_finalize( MPI_COMM_WORLD,IERR )
1206 stop 'Error reading angle file.'
1208 else if (extconf) then
1209 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1210 & write (iout,'(a)') 'Extended chain initial geometry.'
1212 theta(i)=90d0*deg2rad
1215 phi(i)=180d0*deg2rad
1218 alph(i)=110d0*deg2rad
1221 omeg(i)=-120d0*deg2rad
1222 if (itype(i).le.0) omeg(i)=-omeg(i)
1224 call chainbuild_extconf
1226 if(me.eq.king.or..not.out1file)
1227 & write (iout,'(a)') 'Random-generated initial geometry.'
1231 if (me.eq.king .or. fg_rank.eq.0 .and. (
1232 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1236 call gen_rand_conf(itmp,*30)
1238 30 write (iout,*) 'Failed to generate random conformation',
1239 & ', itrial=',itrial
1240 write (*,*) 'Processor:',me,
1241 & ' Failed to generate random conformation',
1251 write (iout,'(a,i3,a)') 'Processor:',me,
1252 & ' error in generating random conformation.'
1253 write (*,'(a,i3,a)') 'Processor:',me,
1254 & ' error in generating random conformation.'
1257 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1262 & ' error in generating random conformation.'
1267 elseif (modecalc.eq.4) then
1268 read (inp,'(a)') intinname
1269 open (intin,file=intinname,status='old',err=333)
1270 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1271 & write (iout,'(a)') 'intinname',intinname
1272 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1274 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1276 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1278 stop 'Error opening angle file.'
1282 C Generate distance constraints, if the PDB structure is to be regularized.
1283 if (nthread.gt.0) then
1284 call read_threadbase
1287 if (me.eq.king .or. .not. out1file)
1289 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1290 write (iout,'(/a,i3,a)')
1291 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1292 write (iout,'(20i4)') (iss(i),i=1,ns)
1294 write(iout,*)"Running with dynamic disulfide-bond formation"
1296 write (iout,'(/a/)') 'Pre-formed links are:'
1302 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1303 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1309 if (ns.gt.0.and.dyn_ss) then
1313 forcon(i-nss)=forcon(i)
1320 dyn_ss_mask(iss(i))=.true.
1325 c write (iout,*) "DC"
1327 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1328 c & (dc(j,i+nres),j=1,3)
1330 if (i2ndstr.gt.0) call secstrp2dihc
1331 c call geom_to_var(nvar,x)
1332 c call etotal(energia(0))
1333 c call enerprint(energia(0))
1334 c call briefout(0,etot)
1336 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1337 cd write (iout,'(a)') 'Variable list:'
1338 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1340 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1341 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1342 & 'Processor',myrank,': end reading molecular data.'
1347 c--------------------------------------------------------------------------
1348 logical function seq_comp(itypea,itypeb,length)
1350 integer length,itypea(length),itypeb(length)
1353 if (itypea(i).ne.itypeb(i)) then
1361 c-----------------------------------------------------------------------------
1362 subroutine read_bridge
1363 C Read information about disulfide bridges.
1364 implicit real*8 (a-h,o-z)
1365 include 'DIMENSIONS'
1369 include 'COMMON.IOUNITS'
1370 include 'COMMON.GEO'
1371 include 'COMMON.VAR'
1372 include 'COMMON.INTERACT'
1373 include 'COMMON.LOCAL'
1374 include 'COMMON.NAMES'
1375 include 'COMMON.CHAIN'
1376 include 'COMMON.FFIELD'
1377 include 'COMMON.SBRIDGE'
1378 include 'COMMON.HEADER'
1379 include 'COMMON.CONTROL'
1380 include 'COMMON.DBASE'
1381 include 'COMMON.THREAD'
1382 include 'COMMON.TIME1'
1383 include 'COMMON.SETUP'
1384 C Read bridging residues.
1385 read (inp,*) ns,(iss(i),i=1,ns)
1387 if(me.eq.king.or..not.out1file)
1388 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1389 C Check whether the specified bridging residues are cystines.
1391 if (itype(iss(i)).ne.1) then
1392 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1393 & 'Do you REALLY think that the residue ',
1394 & restyp(itype(iss(i))),i,
1395 & ' can form a disulfide bridge?!!!'
1396 write (*,'(2a,i3,a)')
1397 & 'Do you REALLY think that the residue ',
1398 & restyp(itype(iss(i))),i,
1399 & ' can form a disulfide bridge?!!!'
1401 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1406 C Read preformed bridges.
1408 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1410 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1413 C Check if the residues involved in bridges are in the specified list of
1414 C bridging residues.
1417 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1418 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1419 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1420 & ' contains residues present in other pairs.'
1421 write (*,'(a,i3,a)') 'Disulfide pair',i,
1422 & ' contains residues present in other pairs.'
1424 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1430 if (ihpb(i).eq.iss(j)) goto 10
1432 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1435 if (jhpb(i).eq.iss(j)) goto 20
1437 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1443 ihpb(i)=ihpb(i)+nres
1444 jhpb(i)=jhpb(i)+nres
1450 c----------------------------------------------------------------------------
1451 subroutine read_x(kanal,*)
1452 implicit real*8 (a-h,o-z)
1453 include 'DIMENSIONS'
1454 include 'COMMON.GEO'
1455 include 'COMMON.VAR'
1456 include 'COMMON.CHAIN'
1457 include 'COMMON.IOUNITS'
1458 include 'COMMON.CONTROL'
1459 include 'COMMON.LOCAL'
1460 include 'COMMON.INTERACT'
1461 c Read coordinates from input
1463 read(kanal,'(8f10.5)',end=10,err=10)
1464 & ((c(l,k),l=1,3),k=1,nres),
1465 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1468 c(j,2*nres)=c(j,nres)
1470 call int_from_cart1(.false.)
1473 dc(j,i)=c(j,i+1)-c(j,i)
1474 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1478 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1480 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1481 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1489 c----------------------------------------------------------------------------
1490 subroutine read_threadbase
1491 implicit real*8 (a-h,o-z)
1492 include 'DIMENSIONS'
1493 include 'COMMON.IOUNITS'
1494 include 'COMMON.GEO'
1495 include 'COMMON.VAR'
1496 include 'COMMON.INTERACT'
1497 include 'COMMON.LOCAL'
1498 include 'COMMON.NAMES'
1499 include 'COMMON.CHAIN'
1500 include 'COMMON.FFIELD'
1501 include 'COMMON.SBRIDGE'
1502 include 'COMMON.HEADER'
1503 include 'COMMON.CONTROL'
1504 include 'COMMON.DBASE'
1505 include 'COMMON.THREAD'
1506 include 'COMMON.TIME1'
1507 C Read pattern database for threading.
1508 read (icbase,*) nseq
1510 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1511 & nres_base(2,i),nres_base(3,i)
1512 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1514 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1515 c & nres_base(2,i),nres_base(3,i)
1516 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1520 if (weidis.eq.0.0D0) weidis=0.1D0
1529 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1530 write (iout,'(a,i5)') 'nexcl: ',nexcl
1531 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1534 c------------------------------------------------------------------------------
1535 subroutine setup_var
1536 implicit real*8 (a-h,o-z)
1537 include 'DIMENSIONS'
1538 include 'COMMON.IOUNITS'
1539 include 'COMMON.GEO'
1540 include 'COMMON.VAR'
1541 include 'COMMON.INTERACT'
1542 include 'COMMON.LOCAL'
1543 include 'COMMON.NAMES'
1544 include 'COMMON.CHAIN'
1545 include 'COMMON.FFIELD'
1546 include 'COMMON.SBRIDGE'
1547 include 'COMMON.HEADER'
1548 include 'COMMON.CONTROL'
1549 include 'COMMON.DBASE'
1550 include 'COMMON.THREAD'
1551 include 'COMMON.TIME1'
1552 C Set up variable list.
1557 write (iout,*) "SETUP_VAR ialph"
1559 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1561 ialph(i,1)=nvar+nside
1565 if (indphi.gt.0) then
1567 else if (indback.gt.0) then
1572 write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1575 c----------------------------------------------------------------------------
1576 subroutine gen_dist_constr
1577 C Generate CA distance constraints.
1578 implicit real*8 (a-h,o-z)
1579 include 'DIMENSIONS'
1580 include 'COMMON.IOUNITS'
1581 include 'COMMON.GEO'
1582 include 'COMMON.VAR'
1583 include 'COMMON.INTERACT'
1584 include 'COMMON.LOCAL'
1585 include 'COMMON.NAMES'
1586 include 'COMMON.CHAIN'
1587 include 'COMMON.FFIELD'
1588 include 'COMMON.SBRIDGE'
1589 include 'COMMON.HEADER'
1590 include 'COMMON.CONTROL'
1591 include 'COMMON.DBASE'
1592 include 'COMMON.THREAD'
1593 include 'COMMON.TIME1'
1594 dimension itype_pdb(maxres)
1595 common /pizda/ itype_pdb
1597 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1598 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1599 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1601 do i=nstart_sup,nstart_sup+nsup-1
1602 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1603 cd & ' seq_pdb', restyp(itype_pdb(i))
1604 do j=i+2,nstart_sup+nsup-1
1606 ihpb(nhpb)=i+nstart_seq-nstart_sup
1607 jhpb(nhpb)=j+nstart_seq-nstart_sup
1609 dhpb(nhpb)=dist(i,j)
1612 cd write (iout,'(a)') 'Distance constraints:'
1617 cd if (ii.gt.nres) then
1622 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1623 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1624 cd & dhpb(i),forcon(i)
1628 c----------------------------------------------------------------------------
1630 implicit real*8 (a-h,o-z)
1631 include 'DIMENSIONS'
1632 include 'COMMON.MAP'
1633 include 'COMMON.IOUNITS'
1634 character*3 angid(4) /'THE','PHI','ALP','OME'/
1635 character*80 mapcard,ucase
1637 read (inp,'(a)') mapcard
1638 mapcard=ucase(mapcard)
1639 if (index(mapcard,'PHI').gt.0) then
1641 else if (index(mapcard,'THE').gt.0) then
1643 else if (index(mapcard,'ALP').gt.0) then
1645 else if (index(mapcard,'OME').gt.0) then
1648 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1649 stop 'Error - illegal variable spec in MAP card.'
1651 call readi (mapcard,'RES1',res1(imap),0)
1652 call readi (mapcard,'RES2',res2(imap),0)
1653 if (res1(imap).eq.0) then
1654 res1(imap)=res2(imap)
1655 else if (res2(imap).eq.0) then
1656 res2(imap)=res1(imap)
1658 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1660 & 'Error - illegal definition of variable group in MAP.'
1661 stop 'Error - illegal definition of variable group in MAP.'
1663 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1664 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1665 call readi(mapcard,'NSTEP',nstep(imap),0)
1666 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1668 & 'Illegal boundary and/or step size specification in MAP.'
1669 stop 'Illegal boundary and/or step size specification in MAP.'
1674 c----------------------------------------------------------------------------
1676 implicit real*8 (a-h,o-z)
1677 include 'DIMENSIONS'
1678 include 'COMMON.IOUNITS'
1679 include 'COMMON.GEO'
1680 include 'COMMON.CSA'
1681 include 'COMMON.BANK'
1682 include 'COMMON.CONTROL'
1684 character*620 mcmcard
1685 call card_concat(mcmcard)
1687 call readi(mcmcard,'NCONF',nconf,50)
1688 call readi(mcmcard,'NADD',nadd,0)
1689 call readi(mcmcard,'JSTART',jstart,1)
1690 call readi(mcmcard,'JEND',jend,1)
1691 call readi(mcmcard,'NSTMAX',nstmax,500000)
1692 call readi(mcmcard,'N0',n0,1)
1693 call readi(mcmcard,'N1',n1,6)
1694 call readi(mcmcard,'N2',n2,4)
1695 call readi(mcmcard,'N3',n3,0)
1696 call readi(mcmcard,'N4',n4,0)
1697 call readi(mcmcard,'N5',n5,0)
1698 call readi(mcmcard,'N6',n6,10)
1699 call readi(mcmcard,'N7',n7,0)
1700 call readi(mcmcard,'N8',n8,0)
1701 call readi(mcmcard,'N9',n9,0)
1702 call readi(mcmcard,'N14',n14,0)
1703 call readi(mcmcard,'N15',n15,0)
1704 call readi(mcmcard,'N16',n16,0)
1705 call readi(mcmcard,'N17',n17,0)
1706 call readi(mcmcard,'N18',n18,0)
1708 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1710 call readi(mcmcard,'NDIFF',ndiff,2)
1711 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1712 call readi(mcmcard,'IS1',is1,1)
1713 call readi(mcmcard,'IS2',is2,8)
1714 call readi(mcmcard,'NRAN0',nran0,4)
1715 call readi(mcmcard,'NRAN1',nran1,2)
1716 call readi(mcmcard,'IRR',irr,1)
1717 call readi(mcmcard,'NSEED',nseed,20)
1718 call readi(mcmcard,'NTOTAL',ntotal,10000)
1719 call reada(mcmcard,'CUT1',cut1,2.0d0)
1720 call reada(mcmcard,'CUT2',cut2,5.0d0)
1721 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1722 call readi(mcmcard,'ICMAX',icmax,3)
1723 call readi(mcmcard,'IRESTART',irestart,0)
1724 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1727 call reada(mcmcard,'DELE',dele,20.0d0)
1728 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1729 call readi(mcmcard,'IREF',iref,0)
1730 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1731 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1732 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1733 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1734 write (iout,*) "NCONF_IN",nconf_in
1737 c----------------------------------------------------------------------------
1738 cfmc subroutine mcmfread
1739 cfmc implicit real*8 (a-h,o-z)
1740 cfmc include 'DIMENSIONS'
1741 cfmc include 'COMMON.MCMF'
1742 cfmc include 'COMMON.IOUNITS'
1743 cfmc include 'COMMON.GEO'
1744 cfmc character*80 ucase
1745 cfmc character*620 mcmcard
1746 cfmc call card_concat(mcmcard)
1748 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1749 cfmc write(iout,*)'MAXRANT=',maxrant
1750 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1751 cfmc write(iout,*)'MAXFAM=',maxfam
1752 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1753 cfmc write(iout,*)'NNET1=',nnet1
1754 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1755 cfmc write(iout,*)'NNET2=',nnet2
1756 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1757 cfmc write(iout,*)'NNET3=',nnet3
1758 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1759 cfmc write(iout,*)'ILASTT=',ilastt
1760 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1761 cfmc write(iout,*)'MAXSTR=',maxstr
1762 cfmc maxstr_f=maxstr/maxfam
1763 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1764 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1765 cfmc write(iout,*)'NMCMF=',nmcmf
1766 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1767 cfmc write(iout,*)'IFOCUS=',ifocus
1768 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1769 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1770 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1771 cfmc write(iout,*)'INTPRT=',intprt
1772 cfmc call readi(mcmcard,'IPRT',iprt,100)
1773 cfmc write(iout,*)'IPRT=',iprt
1774 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1775 cfmc write(iout,*)'IMAXTR=',imaxtr
1776 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1777 cfmc write(iout,*)'MAXEVEN=',maxeven
1778 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1779 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1780 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1781 cfmc write(iout,*)'INIMIN=',inimin
1782 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1783 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1784 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1785 cfmc write(iout,*)'NTHREAD=',nthread
1786 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1787 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1788 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1789 cfmc write(iout,*)'MAXPERT=',maxpert
1790 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1791 cfmc write(iout,*)'IRMSD=',irmsd
1792 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1793 cfmc write(iout,*)'DENEMIN=',denemin
1794 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1795 cfmc write(iout,*)'RCUT1S=',rcut1s
1796 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1797 cfmc write(iout,*)'RCUT1E=',rcut1e
1798 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1799 cfmc write(iout,*)'RCUT2S=',rcut2s
1800 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1801 cfmc write(iout,*)'RCUT2E=',rcut2e
1802 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1803 cfmc write(iout,*)'DPERT1=',d_pert1
1804 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1805 cfmc write(iout,*)'DPERT1A=',d_pert1a
1806 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1807 cfmc write(iout,*)'DPERT2=',d_pert2
1808 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1809 cfmc write(iout,*)'DPERT2A=',d_pert2a
1810 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1811 cfmc write(iout,*)'DPERT2B=',d_pert2b
1812 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1813 cfmc write(iout,*)'DPERT2C=',d_pert2c
1814 cfmc d_pert1=deg2rad*d_pert1
1815 cfmc d_pert1a=deg2rad*d_pert1a
1816 cfmc d_pert2=deg2rad*d_pert2
1817 cfmc d_pert2a=deg2rad*d_pert2a
1818 cfmc d_pert2b=deg2rad*d_pert2b
1819 cfmc d_pert2c=deg2rad*d_pert2c
1820 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1821 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1822 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1823 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1824 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1825 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1826 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1827 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1828 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1829 cfmc write(iout,*)'RCUTINI=',rcutini
1830 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1831 cfmc write(iout,*)'GRAT=',grat
1832 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1833 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1837 c----------------------------------------------------------------------------
1839 implicit real*8 (a-h,o-z)
1840 include 'DIMENSIONS'
1841 include 'COMMON.MCM'
1842 include 'COMMON.MCE'
1843 include 'COMMON.IOUNITS'
1845 character*320 mcmcard
1846 call card_concat(mcmcard)
1847 call readi(mcmcard,'MAXACC',maxacc,100)
1848 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1849 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1850 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1851 call readi(mcmcard,'MAXREPM',maxrepm,200)
1852 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1853 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1854 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1855 call reada(mcmcard,'E_UP',e_up,5.0D0)
1856 call reada(mcmcard,'DELTE',delte,0.1D0)
1857 call readi(mcmcard,'NSWEEP',nsweep,5)
1858 call readi(mcmcard,'NSTEPH',nsteph,0)
1859 call readi(mcmcard,'NSTEPC',nstepc,0)
1860 call reada(mcmcard,'TMIN',tmin,298.0D0)
1861 call reada(mcmcard,'TMAX',tmax,298.0D0)
1862 call readi(mcmcard,'NWINDOW',nwindow,0)
1863 call readi(mcmcard,'PRINT_MC',print_mc,0)
1864 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1865 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1866 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1867 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1868 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1869 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1870 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1871 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1872 if (nwindow.gt.0) then
1873 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1875 winlen(i)=winend(i)-winstart(i)+1
1878 if (tmax.lt.tmin) tmax=tmin
1879 if (tmax.eq.tmin) then
1883 if (nstepc.gt.0 .and. nsteph.gt.0) then
1884 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1885 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1887 C Probabilities of different move types
1888 sumpro_type(0)=0.0D0
1889 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1890 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1891 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1892 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1893 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1894 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1895 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1897 print *,'i',i,' sumprotype',sumpro_type(i)
1898 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1899 print *,'i',i,' sumprotype',sumpro_type(i)
1903 c----------------------------------------------------------------------------
1904 subroutine read_minim
1905 implicit real*8 (a-h,o-z)
1906 include 'DIMENSIONS'
1907 include 'COMMON.MINIM'
1908 include 'COMMON.IOUNITS'
1910 character*320 minimcard
1911 call card_concat(minimcard)
1912 call readi(minimcard,'MAXMIN',maxmin,2000)
1913 call readi(minimcard,'MAXFUN',maxfun,5000)
1914 call readi(minimcard,'MINMIN',minmin,maxmin)
1915 call readi(minimcard,'MINFUN',minfun,maxmin)
1916 call reada(minimcard,'TOLF',tolf,1.0D-2)
1917 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1918 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1919 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1920 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1921 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1922 & 'Options in energy minimization:'
1923 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1924 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1925 & 'MinMin:',MinMin,' MinFun:',MinFun,
1926 & ' TolF:',TolF,' RTolF:',RTolF
1929 c----------------------------------------------------------------------------
1930 subroutine read_angles(kanal,*)
1931 implicit real*8 (a-h,o-z)
1932 include 'DIMENSIONS'
1933 include 'COMMON.GEO'
1934 include 'COMMON.VAR'
1935 include 'COMMON.CHAIN'
1936 include 'COMMON.IOUNITS'
1937 include 'COMMON.CONTROL'
1938 c Read angles from input
1940 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1941 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1942 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1943 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1946 c 9/7/01 avoid 180 deg valence angle
1947 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1949 theta(i)=deg2rad*theta(i)
1950 phi(i)=deg2rad*phi(i)
1951 alph(i)=deg2rad*alph(i)
1952 omeg(i)=deg2rad*omeg(i)
1957 c----------------------------------------------------------------------------
1958 subroutine reada(rekord,lancuch,wartosc,default)
1960 character*(*) rekord,lancuch
1961 double precision wartosc,default
1964 iread=index(rekord,lancuch)
1965 if (iread.eq.0) then
1969 iread=iread+ilen(lancuch)+1
1970 read (rekord(iread:),*,err=10,end=10) wartosc
1975 c----------------------------------------------------------------------------
1976 subroutine readi(rekord,lancuch,wartosc,default)
1978 character*(*) rekord,lancuch
1979 integer wartosc,default
1982 iread=index(rekord,lancuch)
1983 if (iread.eq.0) then
1987 iread=iread+ilen(lancuch)+1
1988 read (rekord(iread:),*,err=10,end=10) wartosc
1993 c----------------------------------------------------------------------------
1994 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1997 integer tablica(dim),default
1998 character*(*) rekord,lancuch
2005 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2006 if (iread.eq.0) return
2007 iread=iread+ilen(lancuch)+1
2008 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2011 c----------------------------------------------------------------------------
2012 subroutine multreada(rekord,lancuch,tablica,dim,default)
2015 double precision tablica(dim),default
2016 character*(*) rekord,lancuch
2023 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2024 if (iread.eq.0) return
2025 iread=iread+ilen(lancuch)+1
2026 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2029 c----------------------------------------------------------------------------
2030 subroutine openunits
2031 implicit real*8 (a-h,o-z)
2032 include 'DIMENSIONS'
2035 character*16 form,nodename
2038 include 'COMMON.SETUP'
2039 include 'COMMON.IOUNITS'
2041 include 'COMMON.CONTROL'
2042 integer lenpre,lenpot,ilen,lentmp
2044 character*3 out1file_text,ucase
2049 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2050 call getenv_loc("PREFIX",prefix)
2052 call getenv_loc("POT",pot)
2053 call getenv_loc("DIRTMP",tmpdir)
2054 call getenv_loc("CURDIR",curdir)
2055 call getenv_loc("OUT1FILE",out1file_text)
2056 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2057 out1file_text=ucase(out1file_text)
2058 if (out1file_text(1:1).eq."Y") then
2061 out1file=fg_rank.gt.0
2066 if (lentmp.gt.0) then
2067 write (*,'(80(1h!))')
2068 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2069 write (*,'(80(1h!))')
2070 write (*,*)"All output files will be on node /tmp directory."
2072 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2073 if (me.eq.king) then
2074 write (*,*) "The master node is ",nodename
2075 else if (fg_rank.eq.0) then
2076 write (*,*) "I am the CG slave node ",nodename
2078 write (*,*) "I am the FG slave node ",nodename
2081 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2082 lenpre = lentmp+lenpre+1
2084 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2085 C Get the names and open the input files
2086 #if defined(WINIFL) || defined(WINPGI)
2087 open(1,file=pref_orig(:ilen(pref_orig))//
2088 & '.inp',status='old',readonly,shared)
2089 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2090 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2091 C Get parameter filenames and open the parameter files.
2092 call getenv_loc('BONDPAR',bondname)
2093 open (ibond,file=bondname,status='old',readonly,shared)
2094 call getenv_loc('THETPAR',thetname)
2095 open (ithep,file=thetname,status='old',readonly,shared)
2096 call getenv_loc('ROTPAR',rotname)
2097 open (irotam,file=rotname,status='old',readonly,shared)
2098 call getenv_loc('TORPAR',torname)
2099 open (itorp,file=torname,status='old',readonly,shared)
2100 call getenv_loc('TORDPAR',tordname)
2101 open (itordp,file=tordname,status='old',readonly,shared)
2102 call getenv_loc('FOURIER',fouriername)
2103 open (ifourier,file=fouriername,status='old',readonly,shared)
2104 call getenv_loc('ELEPAR',elename)
2105 open (ielep,file=elename,status='old',readonly,shared)
2106 call getenv_loc('SIDEPAR',sidename)
2107 open (isidep,file=sidename,status='old',readonly,shared)
2108 call getenv_loc('LIPTRANPAR',liptranname)
2109 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2110 #elif (defined CRAY) || (defined AIX)
2111 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2113 c print *,"Processor",myrank," opened file 1"
2114 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2115 c print *,"Processor",myrank," opened file 9"
2116 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2117 C Get parameter filenames and open the parameter files.
2118 call getenv_loc('BONDPAR',bondname)
2119 open (ibond,file=bondname,status='old',action='read')
2120 c print *,"Processor",myrank," opened file IBOND"
2121 call getenv_loc('THETPAR',thetname)
2122 open (ithep,file=thetname,status='old',action='read')
2124 call getenv_loc('THETPARPDB',thetname_pdb)
2125 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2127 c print *,"Processor",myrank," opened file ITHEP"
2128 call getenv_loc('ROTPAR',rotname)
2129 open (irotam,file=rotname,status='old',action='read')
2131 call getenv_loc('ROTPARPDB',rotname_pdb)
2132 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2134 c print *,"Processor",myrank," opened file IROTAM"
2135 call getenv_loc('TORPAR',torname)
2136 open (itorp,file=torname,status='old',action='read')
2137 c print *,"Processor",myrank," opened file ITORP"
2138 call getenv_loc('TORDPAR',tordname)
2139 open (itordp,file=tordname,status='old',action='read')
2140 c print *,"Processor",myrank," opened file ITORDP"
2141 call getenv_loc('SCCORPAR',sccorname)
2142 open (isccor,file=sccorname,status='old',action='read')
2143 c print *,"Processor",myrank," opened file ISCCOR"
2144 call getenv_loc('FOURIER',fouriername)
2145 open (ifourier,file=fouriername,status='old',action='read')
2146 c print *,"Processor",myrank," opened file IFOURIER"
2147 call getenv_loc('ELEPAR',elename)
2148 open (ielep,file=elename,status='old',action='read')
2149 c print *,"Processor",myrank," opened file IELEP"
2150 call getenv_loc('SIDEPAR',sidename)
2151 open (isidep,file=sidename,status='old',action='read')
2152 call getenv_loc('LIPTRANPAR',liptranname)
2153 open (iliptranpar,file=liptranname,status='old',action='read')
2154 c print *,"Processor",myrank," opened file ISIDEP"
2155 c print *,"Processor",myrank," opened parameter files"
2157 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2158 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2159 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2160 C Get parameter filenames and open the parameter files.
2161 call getenv_loc('BONDPAR',bondname)
2162 open (ibond,file=bondname,status='old')
2163 call getenv_loc('THETPAR',thetname)
2164 open (ithep,file=thetname,status='old')
2166 call getenv_loc('THETPARPDB',thetname_pdb)
2167 open (ithep_pdb,file=thetname_pdb,status='old')
2169 call getenv_loc('ROTPAR',rotname)
2170 open (irotam,file=rotname,status='old')
2172 call getenv_loc('ROTPARPDB',rotname_pdb)
2173 open (irotam_pdb,file=rotname_pdb,status='old')
2175 call getenv_loc('TORPAR',torname)
2176 open (itorp,file=torname,status='old')
2177 call getenv_loc('TORDPAR',tordname)
2178 open (itordp,file=tordname,status='old')
2179 call getenv_loc('SCCORPAR',sccorname)
2180 open (isccor,file=sccorname,status='old')
2181 call getenv_loc('FOURIER',fouriername)
2182 open (ifourier,file=fouriername,status='old')
2183 call getenv_loc('ELEPAR',elename)
2184 open (ielep,file=elename,status='old')
2185 call getenv_loc('SIDEPAR',sidename)
2186 open (isidep,file=sidename,status='old')
2187 call getenv_loc('LIPTRANPAR',liptranname)
2188 open (iliptranpar,file=liptranname,status='old')
2190 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2192 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2193 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2194 C Get parameter filenames and open the parameter files.
2195 call getenv_loc('BONDPAR',bondname)
2196 open (ibond,file=bondname,status='old',readonly)
2197 call getenv_loc('THETPAR',thetname)
2198 open (ithep,file=thetname,status='old',readonly)
2199 call getenv_loc('ROTPAR',rotname)
2200 open (irotam,file=rotname,status='old',readonly)
2201 call getenv_loc('TORPAR',torname)
2202 open (itorp,file=torname,status='old',readonly)
2203 call getenv_loc('TORDPAR',tordname)
2204 open (itordp,file=tordname,status='old',readonly)
2205 call getenv_loc('SCCORPAR',sccorname)
2206 open (isccor,file=sccorname,status='old',readonly)
2208 call getenv_loc('THETPARPDB',thetname_pdb)
2209 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2211 call getenv_loc('FOURIER',fouriername)
2212 open (ifourier,file=fouriername,status='old',readonly)
2213 call getenv_loc('ELEPAR',elename)
2214 open (ielep,file=elename,status='old',readonly)
2215 call getenv_loc('SIDEPAR',sidename)
2216 open (isidep,file=sidename,status='old',readonly)
2217 call getenv_loc('LIPTRANPAR',liptranname)
2218 open (iliptranpar,file=liptranname,status='old',action='read')
2220 call getenv_loc('ROTPARPDB',rotname_pdb)
2221 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2225 call getenv_loc('TUBEPAR',tubename)
2226 #if defined(WINIFL) || defined(WINPGI)
2227 open (itube,file=tubename,status='old',readonly,shared)
2228 #elif (defined CRAY) || (defined AIX)
2229 open (itube,file=tubename,status='old',action='read')
2231 open (itube,file=tubename,status='old')
2233 open (itube,file=tubename,status='old',readonly)
2238 C 8/9/01 In the newest version SCp interaction constants are read from a file
2239 C Use -DOLDSCP to use hard-coded constants instead.
2241 call getenv_loc('SCPPAR',scpname)
2242 #if defined(WINIFL) || defined(WINPGI)
2243 open (iscpp,file=scpname,status='old',readonly,shared)
2244 #elif (defined CRAY) || (defined AIX)
2245 open (iscpp,file=scpname,status='old',action='read')
2247 open (iscpp,file=scpname,status='old')
2249 open (iscpp,file=scpname,status='old',readonly)
2252 call getenv_loc('PATTERN',patname)
2253 #if defined(WINIFL) || defined(WINPGI)
2254 open (icbase,file=patname,status='old',readonly,shared)
2255 #elif (defined CRAY) || (defined AIX)
2256 open (icbase,file=patname,status='old',action='read')
2258 open (icbase,file=patname,status='old')
2260 open (icbase,file=patname,status='old',readonly)
2263 C Open output file only for CG processes
2264 c print *,"Processor",myrank," fg_rank",fg_rank
2265 if (fg_rank.eq.0) then
2267 if (nodes.eq.1) then
2270 npos = dlog10(dfloat(nodes-1))+1
2272 if (npos.lt.3) npos=3
2273 write (liczba,'(i1)') npos
2274 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2276 write (liczba,form) me
2277 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2278 & liczba(:ilen(liczba))
2279 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2281 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2283 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2284 & liczba(:ilen(liczba))//'.mol2'
2285 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2286 & liczba(:ilen(liczba))//'.stat'
2288 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2289 & //liczba(:ilen(liczba))//'.stat')
2290 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2293 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2294 & liczba(:ilen(liczba))//'.const'
2299 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2300 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2301 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2302 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2303 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2305 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2307 rest2name=prefix(:ilen(prefix))//'.rst'
2309 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2312 #if defined(AIX) || defined(PGI) || defined(CRAY)
2313 if (me.eq.king .or. .not. out1file) then
2314 open(iout,file=outname,status='unknown')
2316 open(iout,file="/dev/null",status="unknown")
2320 if (fg_rank.gt.0) then
2321 write (liczba,'(i3.3)') myrank/nfgtasks
2322 write (ll,'(bz,i3.3)') fg_rank
2323 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2329 open(igeom,file=intname,status='unknown',position='append')
2330 open(ipdb,file=pdbname,status='unknown')
2331 open(imol2,file=mol2name,status='unknown')
2332 open(istat,file=statname,status='unknown',position='append')
2334 c1out open(iout,file=outname,status='unknown')
2337 if (me.eq.king .or. .not.out1file)
2338 & open(iout,file=outname,status='unknown')
2340 if (fg_rank.gt.0) then
2341 write (liczba,'(i3.3)') myrank/nfgtasks
2342 write (ll,'(bz,i3.3)') fg_rank
2343 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2348 open(igeom,file=intname,status='unknown',access='append')
2349 open(ipdb,file=pdbname,status='unknown')
2350 open(imol2,file=mol2name,status='unknown')
2351 open(istat,file=statname,status='unknown',access='append')
2353 c1out open(iout,file=outname,status='unknown')
2356 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2357 csa_seed=prefix(:lenpre)//'.CSA.seed'
2358 csa_history=prefix(:lenpre)//'.CSA.history'
2359 csa_bank=prefix(:lenpre)//'.CSA.bank'
2360 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2361 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2362 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2363 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2364 csa_int=prefix(:lenpre)//'.int'
2365 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2366 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2367 csa_in=prefix(:lenpre)//'.CSA.in'
2368 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2371 write (iout,'(80(1h-))')
2372 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2373 write (iout,'(80(1h-))')
2374 write (iout,*) "Input file : ",
2375 & pref_orig(:ilen(pref_orig))//'.inp'
2376 write (iout,*) "Output file : ",
2377 & outname(:ilen(outname))
2379 write (iout,*) "Sidechain potential file : ",
2380 & sidename(:ilen(sidename))
2382 write (iout,*) "SCp potential file : ",
2383 & scpname(:ilen(scpname))
2385 write (iout,*) "Electrostatic potential file : ",
2386 & elename(:ilen(elename))
2387 write (iout,*) "Cumulant coefficient file : ",
2388 & fouriername(:ilen(fouriername))
2389 write (iout,*) "Torsional parameter file : ",
2390 & torname(:ilen(torname))
2391 write (iout,*) "Double torsional parameter file : ",
2392 & tordname(:ilen(tordname))
2393 write (iout,*) "SCCOR parameter file : ",
2394 & sccorname(:ilen(sccorname))
2395 write (iout,*) "Bond & inertia constant file : ",
2396 & bondname(:ilen(bondname))
2397 write (iout,*) "Bending-torsion parameter file : ",
2398 & thetname(:ilen(thetname))
2399 write (iout,*) "Rotamer parameter file : ",
2400 & rotname(:ilen(rotname))
2401 write (iout,*) "Thetpdb parameter file : ",
2402 & thetname_pdb(:ilen(thetname_pdb))
2403 write (iout,*) "Threading database : ",
2404 & patname(:ilen(patname))
2406 &write (iout,*)" DIRTMP : ",
2408 write (iout,'(80(1h-))')
2412 c----------------------------------------------------------------------------
2413 subroutine card_concat(card)
2414 implicit real*8 (a-h,o-z)
2415 include 'DIMENSIONS'
2416 include 'COMMON.IOUNITS'
2418 character*80 karta,ucase
2420 read (inp,'(a)') karta
2423 do while (karta(80:80).eq.'&')
2424 card=card(:ilen(card)+1)//karta(:79)
2425 read (inp,'(a)') karta
2428 card=card(:ilen(card)+1)//karta
2431 c----------------------------------------------------------------------------------
2433 implicit real*8 (a-h,o-z)
2434 include 'DIMENSIONS'
2435 include 'COMMON.CHAIN'
2436 include 'COMMON.IOUNITS'
2438 open(irest2,file=rest2name,status='unknown')
2439 read(irest2,*) totT,EK,potE,totE,t_bath
2442 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2445 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2448 read (irest2,*) iset
2453 c---------------------------------------------------------------------------------
2454 subroutine read_fragments
2455 implicit real*8 (a-h,o-z)
2456 include 'DIMENSIONS'
2460 include 'COMMON.SETUP'
2461 include 'COMMON.CHAIN'
2462 include 'COMMON.IOUNITS'
2464 include 'COMMON.CONTROL'
2465 read(inp,*) nset,nfrag,npair,nfrag_back
2466 loc_qlike=(nfrag_back.lt.0)
2467 nfrag_back=iabs(nfrag_back)
2468 if(me.eq.king.or..not.out1file)
2469 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2470 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2472 read(inp,*) mset(iset)
2474 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2476 if(me.eq.king.or..not.out1file)
2477 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2478 & ifrag(2,i,iset), qinfrag(i,iset)
2481 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2483 if(me.eq.king.or..not.out1file)
2484 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2485 & ipair(2,i,iset), qinpair(i,iset)
2489 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2490 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2491 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2492 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2493 if(me.eq.king.or..not.out1file)
2494 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2495 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2496 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2497 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2501 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2502 & wfrag_back(3,i,iset),
2503 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2504 if(me.eq.king.or..not.out1file)
2505 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2506 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2512 C---------------------------------------------------------------------------
2513 subroutine read_afminp
2514 implicit real*8 (a-h,o-z)
2515 include 'DIMENSIONS'
2519 include 'COMMON.SETUP'
2520 include 'COMMON.CONTROL'
2521 include 'COMMON.CHAIN'
2522 include 'COMMON.IOUNITS'
2523 include 'COMMON.SBRIDGE'
2524 character*320 afmcard
2526 call card_concat(afmcard)
2527 call readi(afmcard,"BEG",afmbeg,0)
2528 call readi(afmcard,"END",afmend,0)
2529 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2530 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2531 print *,'FORCE=' ,forceAFMconst
2532 CCCC NOW PROPERTIES FOR AFM
2535 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2537 distafminit=dsqrt(distafminit)
2538 print *,'initdist',distafminit
2541 c-------------------------------------------------------------------------------
2542 subroutine read_saxs_constr
2543 implicit real*8 (a-h,o-z)
2544 include 'DIMENSIONS'
2548 include 'COMMON.SETUP'
2549 include 'COMMON.CONTROL'
2550 include 'COMMON.CHAIN'
2551 include 'COMMON.IOUNITS'
2552 include 'COMMON.SBRIDGE'
2553 double precision cm(3)
2555 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2557 if (saxs_mode.eq.0) then
2558 c SAXS distance distribution
2560 read(inp,*) distsaxs(i),Psaxs(i)
2564 Cnorm = Cnorm + Psaxs(i)
2566 write (iout,*) "Cnorm",Cnorm
2568 Psaxs(i)=Psaxs(i)/Cnorm
2570 write (iout,*) "Normalized distance distribution from SAXS"
2572 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2576 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2578 write (iout,*) "Wsaxs0",Wsaxs0
2582 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2589 cm(j)=cm(j)+Csaxs(j,i)
2597 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2600 write (iout,*) "SAXS sphere coordinates"
2602 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2608 c-------------------------------------------------------------------------------
2609 subroutine read_dist_constr
2610 implicit real*8 (a-h,o-z)
2611 include 'DIMENSIONS'
2615 include 'COMMON.SETUP'
2616 include 'COMMON.CONTROL'
2617 include 'COMMON.CHAIN'
2618 include 'COMMON.IOUNITS'
2619 include 'COMMON.SBRIDGE'
2620 include 'COMMON.INTERACT'
2621 integer ifrag_(2,100),ipair_(2,1000)
2622 double precision wfrag_(100),wpair_(1000)
2623 character*5000 controlcard
2624 logical normalize,next
2626 double precision scal_bfac
2627 double precision xlink(4,0:4) /
2629 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2630 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2631 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2632 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2633 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2634 c print *, "WCHODZE"
2635 write (iout,*) "Calling read_dist_constr"
2636 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2638 restr_on_coord=.false.
2647 call card_concat(controlcard)
2648 next = index(controlcard,"NEXT").gt.0
2649 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2650 write (iout,*) "restr_type",restr_type
2651 call readi(controlcard,"NFRAG",nfrag_,0)
2652 call readi(controlcard,"NFRAG",nfrag_,0)
2653 call readi(controlcard,"NPAIR",npair_,0)
2654 call readi(controlcard,"NDIST",ndist_,0)
2655 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2656 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2657 if (restr_type.eq.10)
2658 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2659 if (restr_type.eq.12)
2660 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2661 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2662 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2663 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2664 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2665 normalize = index(controlcard,"NORMALIZE").gt.0
2666 write (iout,*) "WBOLTZD",wboltzd
2667 write (iout,*) "SCAL_PEAK",scal_peak
2668 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2669 write (iout,*) "IFRAG"
2671 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2673 write (iout,*) "IPAIR"
2675 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2677 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2679 & "Distance restraints as generated from reference structure"
2681 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2682 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2683 & ifrag_(2,i)=nstart_sup+nsup-1
2684 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2686 if (wfrag_(i).eq.0.0d0) cycle
2687 do j=ifrag_(1,i),ifrag_(2,i)-1
2688 do k=j+1,ifrag_(2,i)
2689 c write (iout,*) "j",j," k",k
2691 if (restr_type.eq.1) then
2697 forcon(nhpb)=wfrag_(i)
2698 else if (constr_dist.eq.2) then
2699 if (ddjk.le.dist_cut) then
2705 forcon(nhpb)=wfrag_(i)
2707 else if (restr_type.eq.3) then
2713 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2716 if (.not.out1file .or. me.eq.king)
2717 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2718 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2720 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2721 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2727 if (wpair_(i).eq.0.0d0) cycle
2735 do j=ifrag_(1,ii),ifrag_(2,ii)
2736 do k=ifrag_(1,jj),ifrag_(2,jj)
2738 if (restr_type.eq.1) then
2744 forcon(nhpb)=wpair_(i)
2745 else if (constr_dist.eq.2) then
2746 if (ddjk.le.dist_cut) then
2752 forcon(nhpb)=wpair_(i)
2754 else if (restr_type.eq.3) then
2760 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2763 if (.not.out1file .or. me.eq.king)
2764 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2765 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2767 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2768 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2775 write (iout,*) "Distance restraints as read from input"
2777 if (restr_type.eq.12) then
2778 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2779 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2780 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2781 & fordepth_peak(nhpb_peak+1),npeak
2782 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2783 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2784 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2785 c & fordepth_peak(nhpb_peak+1),npeak
2786 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2787 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2788 nhpb_peak=nhpb_peak+1
2789 irestr_type_peak(nhpb_peak)=12
2790 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2793 if (.not.out1file .or. me.eq.king)
2794 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2795 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2796 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2797 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2798 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2800 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2801 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2802 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2803 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2804 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2806 if (ibecarb_peak(nhpb_peak).eq.3) then
2807 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2808 else if (ibecarb_peak(nhpb_peak).eq.2) then
2809 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2810 else if (ibecarb_peak(nhpb_peak).eq.1) then
2811 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2812 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2814 else if (restr_type.eq.11) then
2815 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2816 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2817 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2818 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2820 irestr_type(nhpb)=11
2822 if (.not.out1file .or. me.eq.king)
2823 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2824 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2825 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2827 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2828 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2829 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2831 c if (ibecarb(nhpb).gt.0) then
2832 c ihpb(nhpb)=ihpb(nhpb)+nres
2833 c jhpb(nhpb)=jhpb(nhpb)+nres
2835 if (ibecarb(nhpb).eq.3) then
2836 ihpb(nhpb)=ihpb(nhpb)+nres
2837 else if (ibecarb(nhpb).eq.2) then
2838 ihpb(nhpb)=ihpb(nhpb)+nres
2839 else if (ibecarb(nhpb).eq.1) then
2840 ihpb(nhpb)=ihpb(nhpb)+nres
2841 jhpb(nhpb)=jhpb(nhpb)+nres
2843 else if (restr_type.eq.10) then
2844 c Cross-lonk Markov-like potential
2845 call card_concat(controlcard)
2846 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2847 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2849 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2850 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2851 if (index(controlcard,"ZL").gt.0) then
2853 else if (index(controlcard,"ADH").gt.0) then
2855 else if (index(controlcard,"PDH").gt.0) then
2857 else if (index(controlcard,"DSS").gt.0) then
2862 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2863 & xlink(1,link_type))
2864 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2865 & xlink(2,link_type))
2866 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2867 & xlink(3,link_type))
2868 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2869 & xlink(4,link_type))
2870 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2871 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2872 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2873 if (forcon(nhpb+1).le.0.0d0 .or.
2874 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2876 irestr_type(nhpb)=10
2877 if (ibecarb(nhpb).eq.3) then
2878 jhpb(nhpb)=jhpb(nhpb)+nres
2879 else if (ibecarb(nhpb).eq.2) then
2880 ihpb(nhpb)=ihpb(nhpb)+nres
2881 else if (ibecarb(nhpb).eq.1) then
2882 ihpb(nhpb)=ihpb(nhpb)+nres
2883 jhpb(nhpb)=jhpb(nhpb)+nres
2886 if (.not.out1file .or. me.eq.king)
2887 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2888 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2889 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2892 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2893 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2894 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2899 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2900 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2901 if (forcon(nhpb+1).gt.0.0d0) then
2903 if (dhpb1(nhpb).eq.0.0d0) then
2908 if (ibecarb(nhpb).eq.3) then
2909 jhpb(nhpb)=jhpb(nhpb)+nres
2910 else if (ibecarb(nhpb).eq.2) then
2911 ihpb(nhpb)=ihpb(nhpb)+nres
2912 else if (ibecarb(nhpb).eq.1) then
2913 ihpb(nhpb)=ihpb(nhpb)+nres
2914 jhpb(nhpb)=jhpb(nhpb)+nres
2916 if (dhpb(nhpb).eq.0.0d0)
2917 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2920 if (.not.out1file .or. me.eq.king)
2921 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2922 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2924 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2925 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2928 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2929 C if (forcon(nhpb+1).gt.0.0d0) then
2931 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2934 if (restr_type.eq.4) then
2935 write (iout,*) "The BFAC array"
2937 write (iout,'(i5,f10.5)') i,bfac(i)
2940 if (itype(i).eq.ntyp1) cycle
2942 if (itype(j).eq.ntyp1) cycle
2943 if (itype(i).eq.10) then
2948 if (itype(j).eq.10) then
2958 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2961 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2962 ihpb(nhpb)=i+nres*ii
2963 jhpb(nhpb)=j+nres*jj
2964 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2966 if (.not.out1file .or. me.eq.king) then
2967 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2968 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2969 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2973 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2974 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2975 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2985 if (restr_type.eq.5) then
2986 restr_on_coord=.true.
2988 if (itype(i).eq.ntyp1) cycle
2989 bfac(i)=(scal_bfac/bfac(i))**2
2998 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2999 & fordepthmax=fordepth(i)
3002 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3007 c-------------------------------------------------------------------------------
3008 subroutine read_constr_homology
3010 include 'DIMENSIONS'
3014 include 'COMMON.SETUP'
3015 include 'COMMON.CONTROL'
3016 include 'COMMON.CHAIN'
3017 include 'COMMON.IOUNITS'
3019 include 'COMMON.GEO'
3020 include 'COMMON.INTERACT'
3021 include 'COMMON.NAMES'
3023 c For new homol impl
3025 include 'COMMON.VAR'
3028 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3030 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3031 c & sigma_odl_temp(maxres,maxres,max_template)
3033 character*24 model_ki_dist, model_ki_angle
3034 character*500 controlcard
3035 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3040 c FP - Nov. 2014 Temporary specifications for new vars
3042 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3044 double precision, dimension (max_template,maxres) :: rescore
3045 double precision, dimension (max_template,maxres) :: rescore2
3046 double precision, dimension (max_template,maxres) :: rescore3
3047 character*24 pdbfile,tpl_k_rescore
3048 c -----------------------------------------------------------------
3049 c Reading multiple PDB ref structures and calculation of retraints
3050 c not using pre-computed ones stored in files model_ki_{dist,angle}
3052 c -----------------------------------------------------------------
3055 c Alternative: reading from input
3056 call card_concat(controlcard)
3057 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3058 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3059 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3060 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3061 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3062 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3063 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3064 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3065 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3066 if(.not.read2sigma.and.start_from_model) then
3067 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3068 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3069 start_from_model=.false.
3071 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3072 & write(iout,*) 'START_FROM_MODELS is ON'
3073 if(start_from_model .and. rest) then
3074 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3075 write(iout,*) 'START_FROM_MODELS is OFF'
3076 write(iout,*) 'remove restart keyword from input'
3079 if (homol_nset.gt.1)then
3080 call card_concat(controlcard)
3081 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3082 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3083 write(iout,*) "iset homology_weight "
3085 write(iout,*) i,waga_homology(i)
3088 iset=mod(kolor,homol_nset)+1
3091 waga_homology(1)=1.0
3094 cd write (iout,*) "nnt",nnt," nct",nct
3101 c write(iout,*) 'nnt=',nnt,'nct=',nct
3104 do k=1,constr_homology
3117 if (read_homol_frag) then
3118 call read_klapaucjusz
3121 do k=1,constr_homology
3123 read(inp,'(a)') pdbfile
3124 if(me.eq.king .or. .not. out1file)
3125 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3126 & pdbfile(:ilen(pdbfile))
3127 open(ipdbin,file=pdbfile,status='old',err=33)
3129 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3130 & pdbfile(:ilen(pdbfile))
3133 c print *,'Begin reading pdb data'
3135 c Files containing res sim or local scores (former containing sigmas)
3138 write(kic2,'(bz,i2.2)') k
3140 tpl_k_rescore="template"//kic2//".sco"
3143 if (read2sigma) then
3144 call readpdb_template(k)
3149 c Distance restraints
3152 C Copy the coordinates from reference coordinates (?)
3156 c write (iout,*) "c(",j,i,") =",c(j,i)
3160 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3162 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3163 open (ientin,file=tpl_k_rescore,status='old')
3164 if (nnt.gt.1) rescore(k,1)=0.0d0
3165 do irec=nnt,nct ! loop for reading res sim
3166 if (read2sigma) then
3167 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3168 & rescore3_tmp,idomain_tmp
3170 idomain(k,i_tmp)=idomain_tmp
3171 rescore(k,i_tmp)=rescore_tmp
3172 rescore2(k,i_tmp)=rescore2_tmp
3173 rescore3(k,i_tmp)=rescore3_tmp
3174 if (.not. out1file .or. me.eq.king)
3175 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3176 & i_tmp,rescore2_tmp,rescore_tmp,
3177 & rescore3_tmp,idomain_tmp
3180 read (ientin,*,end=1401) rescore_tmp
3182 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3183 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3184 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3189 if (waga_dist.ne.0.0d0) then
3197 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3198 c write (iout,*) k,i,j,distal,dist2_cut
3200 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3201 & .and. distal.le.dist2_cut ) then
3207 c write (iout,*) "k",k
3208 c write (iout,*) "i",i," j",j," constr_homology",
3213 if (read2sigma) then
3216 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3218 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3219 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3220 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3222 if (odl(k,ii).le.dist_cut) then
3223 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3226 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3227 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3229 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3230 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3234 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3237 l_homo(k,ii)=.false.
3244 c Theta, dihedral and SC retraints
3246 if (waga_angle.gt.0.0d0) then
3247 c open (ientin,file=tpl_k_sigma_dih,status='old')
3248 c do irec=1,maxres-3 ! loop for reading sigma_dih
3249 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3250 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3251 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3252 c & sigma_dih(k,i+nnt-1)
3257 if (idomain(k,i).eq.0) then
3261 dih(k,i)=phiref(i) ! right?
3262 c read (ientin,*) sigma_dih(k,i) ! original variant
3263 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3264 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3265 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3266 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3267 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3269 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3270 & rescore(k,i-2)+rescore(k,i-3))/4.0
3271 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3272 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3273 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3274 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3275 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3276 c Instead of res sim other local measure of b/b str reliability possible
3277 if (sigma_dih(k,i).ne.0)
3278 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3279 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3284 if (waga_theta.gt.0.0d0) then
3285 c open (ientin,file=tpl_k_sigma_theta,status='old')
3286 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3287 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3288 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3289 c & sigma_theta(k,i+nnt-1)
3294 do i = nnt+2,nct ! right? without parallel.
3295 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3296 c do i=ithet_start,ithet_end ! with FG parallel.
3297 if (idomain(k,i).eq.0) then
3298 sigma_theta(k,i)=0.0
3301 thetatpl(k,i)=thetaref(i)
3302 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3303 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3304 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3305 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3306 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3307 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3308 & rescore(k,i-2))/3.0
3309 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3310 if (sigma_theta(k,i).ne.0)
3311 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3313 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3314 c rescore(k,i-2) ! right expression ?
3315 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3319 if (waga_d.gt.0.0d0) then
3320 c open (ientin,file=tpl_k_sigma_d,status='old')
3321 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3322 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3323 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3324 c & sigma_d(k,i+nnt-1)
3328 do i = nnt,nct ! right? without parallel.
3329 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3330 c do i=loc_start,loc_end ! with FG parallel.
3331 if (itype(i).eq.10) cycle
3332 if (idomain(k,i).eq.0 ) then
3339 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3340 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3341 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3342 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3343 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3344 if (sigma_d(k,i).ne.0)
3345 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3347 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3348 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3349 c read (ientin,*) sigma_d(k,i) ! 1st variant
3354 c remove distance restraints not used in any model from the list
3355 c shift data in all arrays
3357 if (waga_dist.ne.0.0d0) then
3363 if (ii_in_use(ii).eq.0.and.liiflag) then
3367 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3368 & .not.liiflag.and.ii.eq.lim_odl) then
3369 if (ii.eq.lim_odl) then
3370 iishift=ii-iistart+1
3375 do ki=iistart,lim_odl-iishift
3376 ires_homo(ki)=ires_homo(ki+iishift)
3377 jres_homo(ki)=jres_homo(ki+iishift)
3378 ii_in_use(ki)=ii_in_use(ki+iishift)
3379 do k=1,constr_homology
3380 odl(k,ki)=odl(k,ki+iishift)
3381 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3382 l_homo(k,ki)=l_homo(k,ki+iishift)
3386 lim_odl=lim_odl-iishift
3392 endif ! .not. klapaucjusz
3394 if (constr_homology.gt.0) call homology_partition
3395 if (constr_homology.gt.0) call init_int_table
3396 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3397 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3401 if (.not.out_template_restr) return
3402 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3403 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3404 write (iout,*) "Distance restraints from templates"
3406 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3407 & ii,ires_homo(ii),jres_homo(ii),
3408 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3409 & ki=1,constr_homology)
3411 write (iout,*) "Dihedral angle restraints from templates"
3413 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3414 & (rad2deg*dih(ki,i),
3415 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3417 write (iout,*) "Virtual-bond angle restraints from templates"
3419 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3420 & (rad2deg*thetatpl(ki,i),
3421 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3423 write (iout,*) "SC restraints from templates"
3425 write(iout,'(i5,100(4f8.2,4x))') i,
3426 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3427 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3430 c -----------------------------------------------------------------
3433 c----------------------------------------------------------------------
3435 subroutine flush(iu)
3440 subroutine flush(iu)
3445 c------------------------------------------------------------------------------
3446 subroutine copy_to_tmp(source)
3447 include "DIMENSIONS"
3448 include "COMMON.IOUNITS"
3449 character*(*) source
3450 character* 256 tmpfile
3454 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3455 inquire(file=tmpfile,exist=ex)
3457 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3458 & " to temporary directory..."
3459 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3460 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3464 c------------------------------------------------------------------------------
3465 subroutine move_from_tmp(source)
3466 include "DIMENSIONS"
3467 include "COMMON.IOUNITS"
3468 character*(*) source
3471 write (*,*) "Moving ",source(:ilen(source)),
3472 & " from temporary directory to working directory"
3473 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3474 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3477 c------------------------------------------------------------------------------
3478 subroutine random_init(seed)
3480 C Initialize random number generator
3482 implicit real*8 (a-h,o-z)
3483 include 'DIMENSIONS'
3486 logical OKRandom, prng_restart
3488 integer iseed_array(4)
3490 include 'COMMON.IOUNITS'
3491 include 'COMMON.TIME1'
3492 include 'COMMON.THREAD'
3493 include 'COMMON.SBRIDGE'
3494 include 'COMMON.CONTROL'
3495 include 'COMMON.MCM'
3496 include 'COMMON.MAP'
3497 include 'COMMON.HEADER'
3498 include 'COMMON.CSA'
3499 include 'COMMON.CHAIN'
3500 include 'COMMON.MUCA'
3502 include 'COMMON.FFIELD'
3503 include 'COMMON.SETUP'
3504 iseed=-dint(dabs(seed))
3505 if (iseed.eq.0) then
3506 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3507 & 'Random seed undefined. The program will stop.'
3508 write (*,'(/80(1h*)/20x,a/80(1h*))')
3509 & 'Random seed undefined. The program will stop.'
3511 call mpi_finalize(mpi_comm_world,ierr)
3513 stop 'Bad random seed.'
3516 if (fg_rank.eq.0) then
3520 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3521 OKRandom = prng_restart(me,iseed)
3524 tmp=65536.0d0**(4-i)
3525 iseed_array(i) = dint(seed/tmp)
3526 seed=seed-iseed_array(i)*tmp
3529 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3530 & (iseed_array(i),i=1,4)
3531 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3532 & (iseed_array(i),i=1,4)
3533 OKRandom = prng_restart(me,iseed_array)
3536 c r1 = prng_next(me)
3537 r1=ran_number(0.0D0,1.0D0)
3539 & write (iout,*) 'ran_num',r1
3540 if (r1.lt.0.0d0) OKRandom=.false.
3542 if (.not.OKRandom) then
3543 write (iout,*) 'PRNG IS NOT WORKING!!!'
3544 print *,'PRNG IS NOT WORKING!!!'
3547 call mpi_abort(mpi_comm_world,error_msg,ierr)
3550 write (iout,*) 'too many processors for parallel prng'
3551 write (*,*) 'too many processors for parallel prng'
3559 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3563 c----------------------------------------------------------------------
3564 subroutine read_klapaucjusz
3566 include 'DIMENSIONS'
3570 include 'COMMON.SETUP'
3571 include 'COMMON.CONTROL'
3572 include 'COMMON.CHAIN'
3573 include 'COMMON.IOUNITS'
3575 include 'COMMON.GEO'
3576 include 'COMMON.INTERACT'
3577 include 'COMMON.NAMES'
3578 character*256 fragfile
3579 integer ninclust(maxclust),inclust(max_template,maxclust),
3580 & nresclust(maxclust),iresclust(maxres,maxclust)
3583 character*24 model_ki_dist, model_ki_angle
3584 character*500 controlcard
3585 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3586 logical lprn /.true./
3592 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3593 double precision, dimension (max_template,maxres) :: rescore
3594 double precision, dimension (max_template,maxres) :: rescore2
3595 character*24 pdbfile,tpl_k_rescore
3598 c For new homol impl
3600 include 'COMMON.VAR'
3602 call getenv("FRAGFILE",fragfile)
3603 open(ientin,file=fragfile,status="old",err=10)
3604 read(ientin,*) constr_homology,nclust
3610 do k=1,constr_homology
3611 read(ientin,'(a)') pdbfile
3612 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3613 & pdbfile(:ilen(pdbfile))
3614 open(ipdbin,file=pdbfile,status='old',err=33)
3616 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3617 & pdbfile(:ilen(pdbfile))
3621 call readpdb_template(k)
3629 read(ientin,*) ninclust(i),nresclust(i)
3630 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3631 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3634 c Loop over clusters
3637 do ll = 1,ninclust(l)
3645 idomain(k,iresclust(i,l)+1) = 1
3647 idomain(k,iresclust(i,l)) = 1
3651 c Distance restraints
3654 C Copy the coordinates from reference coordinates (?)
3658 c write (iout,*) "c(",j,i,") =",c(j,i)
3661 call int_from_cart(.true.,.false.)
3662 call sc_loc_geom(.false.)
3664 thetaref(i)=theta(i)
3667 if (waga_dist.ne.0.0d0) then
3675 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3676 c write (iout,*) k,i,j,distal,dist2_cut
3678 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3679 & .and. distal.le.dist2_cut ) then
3685 c write (iout,*) "k",k
3686 c write (iout,*) "i",i," j",j," constr_homology",
3691 if (read2sigma) then
3694 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3696 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3697 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3698 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3700 if (odl(k,ii).le.dist_cut) then
3701 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3704 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3705 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3707 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3708 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3712 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3715 c l_homo(k,ii)=.false.
3722 c Theta, dihedral and SC retraints
3724 if (waga_angle.gt.0.0d0) then
3726 if (idomain(k,i).eq.0) then
3727 c sigma_dih(k,i)=0.0
3731 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3732 & rescore(k,i-2)+rescore(k,i-3))/4.0
3733 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3734 c & " sigma_dihed",sigma_dih(k,i)
3735 if (sigma_dih(k,i).ne.0)
3736 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3741 if (waga_theta.gt.0.0d0) then
3743 if (idomain(k,i).eq.0) then
3744 c sigma_theta(k,i)=0.0
3747 thetatpl(k,i)=thetaref(i)
3748 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3749 & rescore(k,i-2))/3.0
3750 if (sigma_theta(k,i).ne.0)
3751 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3755 if (waga_d.gt.0.0d0) then
3757 if (itype(i).eq.10) cycle
3758 if (idomain(k,i).eq.0 ) then
3765 sigma_d(k,i)=rescore(k,i)
3766 if (sigma_d(k,i).ne.0)
3767 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3768 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3774 c remove distance restraints not used in any model from the list
3775 c shift data in all arrays
3777 if (waga_dist.ne.0.0d0) then
3783 if (ii_in_use(ii).eq.0.and.liiflag) then
3787 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3788 & .not.liiflag.and.ii.eq.lim_odl) then
3789 if (ii.eq.lim_odl) then
3790 iishift=ii-iistart+1
3795 do ki=iistart,lim_odl-iishift
3796 ires_homo(ki)=ires_homo(ki+iishift)
3797 jres_homo(ki)=jres_homo(ki+iishift)
3798 ii_in_use(ki)=ii_in_use(ki+iishift)
3799 do k=1,constr_homology
3800 odl(k,ki)=odl(k,ki+iishift)
3801 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3802 l_homo(k,ki)=l_homo(k,ki+iishift)
3806 lim_odl=lim_odl-iishift
3813 10 stop "Error infragment file"