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 molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 include 'COMMON.SHIELD'
86 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
87 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
89 character*320 controlcard
94 read (INP,'(a)') titel
95 call card_concat(controlcard)
96 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
97 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
98 call reada(controlcard,'SEED',seed,0.0D0)
99 call random_init(seed)
100 C Set up the time limit (caution! The time must be input in minutes!)
101 read_cart=index(controlcard,'READ_CART').gt.0
102 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
103 C this variable with_theta_constr is the variable which allow to read and execute the
104 C constrains on theta angles WITH_THETA_CONSTR is the keyword
105 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
106 write (iout,*) "with_theta_constr ",with_theta_constr
107 call readi(controlcard,'SYM',symetr,1)
108 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
109 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
110 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
111 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
112 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
113 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
114 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
115 call reada(controlcard,'DRMS',drms,0.1D0)
116 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
117 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
118 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
119 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
120 write (iout,'(a,f10.1)')'DRMS = ',drms
121 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
122 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
124 call readi(controlcard,'NZ_START',nz_start,0)
125 call readi(controlcard,'NZ_END',nz_end,0)
126 call readi(controlcard,'IZ_SC',iz_sc,0)
128 safety = 60.0d0*safety
131 call reada(controlcard,"T_BATH",t_bath,300.0d0)
132 minim=(index(controlcard,'MINIMIZE').gt.0)
133 dccart=(index(controlcard,'CART').gt.0)
134 overlapsc=(index(controlcard,'OVERLAP').gt.0)
135 overlapsc=.not.overlapsc
136 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
137 searchsc=.not.searchsc
138 sideadd=(index(controlcard,'SIDEADD').gt.0)
139 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
140 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
141 outpdb=(index(controlcard,'PDBOUT').gt.0)
142 outmol2=(index(controlcard,'MOL2OUT').gt.0)
143 pdbref=(index(controlcard,'PDBREF').gt.0)
144 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
145 call readi(controlcard,'PDBSTART',indpdb,0)
146 if (index(controlcard,'PDBSTART').gt.0 .and. indpdb.eq.0) indpdb=2
147 write (iout,*) "indpdb",indpdb
148 extconf=(index(controlcard,'EXTCONF').gt.0)
149 AFMlog=(index(controlcard,'AFM'))
150 selfguide=(index(controlcard,'SELFGUIDE'))
151 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
152 call readi(controlcard,'TUBEMOD',tubelog,0)
153 c write (iout,*) TUBElog,"TUBEMODE"
154 call readi(controlcard,'IPRINT',iprint,0)
155 C SHIELD keyword sets if the shielding effect of side-chains is used
156 C 0 denots no shielding is used all peptide are equally despite the
157 C solvent accesible area
158 C 1 the newly introduced function
159 C 2 reseved for further possible developement
160 call readi(controlcard,'SHIELD',shield_mode,0)
161 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
162 write(iout,*) "shield_mode",shield_mode
164 call readi(controlcard,'TORMODE',tor_mode,0)
165 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
166 write(iout,*) "torsional and valence angle mode",tor_mode
167 call readi(controlcard,'MAXGEN',maxgen,10000)
168 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
169 call readi(controlcard,"KDIAG",kdiag,0)
170 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
171 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
172 & write (iout,*) "RESCALE_MODE",rescale_mode
173 split_ene=index(controlcard,'SPLIT_ENE').gt.0
174 if (index(controlcard,'REGULAR').gt.0.0D0) then
175 call reada(controlcard,'WEIDIS',weidis,0.1D0)
179 if (index(controlcard,'CHECKGRAD').gt.0) then
181 if (index(controlcard,' CART').gt.0) then
183 elseif (index(controlcard,'CARINT').gt.0) then
188 call reada(controlcard,'DELTA',aincr,1.0d-7)
189 c write (iout,*) "icheckgrad",icheckgrad
190 elseif (index(controlcard,'THREAD').gt.0) then
192 call readi(controlcard,'THREAD',nthread,0)
193 if (nthread.gt.0) then
194 call reada(controlcard,'WEIDIS',weidis,0.1D0)
197 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
198 stop 'Error termination in Read_Control.'
200 else if (index(controlcard,'MCMA').gt.0) then
202 else if (index(controlcard,'MCEE').gt.0) then
204 else if (index(controlcard,'MULTCONF').gt.0) then
206 else if (index(controlcard,'MAP').gt.0) then
208 call readi(controlcard,'MAP',nmap,0)
209 else if (index(controlcard,'CSA').gt.0) then
211 crc else if (index(controlcard,'ZSCORE').gt.0) then
213 crc ZSCORE is rm from UNRES, modecalc=9 is available
216 cfcm else if (index(controlcard,'MCMF').gt.0) then
218 else if (index(controlcard,'SOFTREG').gt.0) then
220 else if (index(controlcard,'CHECK_BOND').gt.0) then
222 else if (index(controlcard,'TEST').gt.0) then
224 else if (index(controlcard,'MD').gt.0) then
226 else if (index(controlcard,'RE ').gt.0) then
230 lmuca=index(controlcard,'MUCA').gt.0
231 call readi(controlcard,'MUCADYN',mucadyn,0)
232 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
233 if (lmuca .and. (me.eq.king .or. .not.out1file ))
235 write (iout,*) 'MUCADYN=',mucadyn
236 write (iout,*) 'MUCASMOOTH=',muca_smooth
239 iscode=index(controlcard,'ONE_LETTER')
240 indphi=index(controlcard,'PHI')
241 indback=index(controlcard,'BACK')
242 iranconf=index(controlcard,'RAND_CONF')
243 i2ndstr=index(controlcard,'USE_SEC_PRED')
244 gradout=index(controlcard,'GRADOUT').gt.0
245 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
246 C DISTCHAINMAX become obsolete for periodic boundry condition
247 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
248 C Reading the dimensions of box in x,y,z coordinates
249 call reada(controlcard,'BOXX',boxxsize,100.0d0)
250 call reada(controlcard,'BOXY',boxysize,100.0d0)
251 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
252 c Cutoff range for interactions
253 call reada(controlcard,"R_CUT",r_cut,15.0d0)
254 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
255 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
256 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
257 if (lipthick.gt.0.0d0) then
258 bordliptop=(boxzsize+lipthick)/2.0
259 bordlipbot=bordliptop-lipthick
261 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
262 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
263 buflipbot=bordlipbot+lipbufthick
264 bufliptop=bordliptop-lipbufthick
265 if ((lipbufthick*2.0d0).gt.lipthick)
266 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
268 write(iout,*) "bordliptop=",bordliptop
269 write(iout,*) "bordlipbot=",bordlipbot
270 write(iout,*) "bufliptop=",bufliptop
271 write(iout,*) "buflipbot=",buflipbot
272 write (iout,*) "SHIELD MODE",shield_mode
273 if (TUBElog.gt.0) then
274 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
275 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
276 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
277 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
278 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
279 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
280 buftubebot=bordtubebot+tubebufthick
281 buftubetop=bordtubetop-tubebufthick
283 c if (shield_mode.gt.0) then
285 C VSolvSphere the volume of solving sphere
287 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
288 C there will be no distinction between proline peptide group and normal peptide
289 C group in case of shielding parameters
290 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
291 c VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
292 c VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
293 c write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
295 C long axis of side chain
297 c long_r_sidechain(i)=vbldsc0(1,i)
298 c short_r_sidechain(i)=sigma0(i)
302 if (me.eq.king .or. .not.out1file )
303 & write (iout,*) "DISTCHAINMAX",distchainmax
305 if(me.eq.king.or..not.out1file)
306 & write (iout,'(2a)') diagmeth(kdiag),
307 & ' routine used to diagonalize matrices.'
310 c--------------------------------------------------------------------------
311 subroutine read_REMDpar
315 implicit real*8 (a-h,o-z)
317 include 'COMMON.IOUNITS'
318 include 'COMMON.TIME1'
321 include 'COMMON.LANGEVIN'
323 include 'COMMON.LANGEVIN.lang0'
325 include 'COMMON.INTERACT'
326 include 'COMMON.NAMES'
328 include 'COMMON.REMD'
329 include 'COMMON.CONTROL'
330 include 'COMMON.SETUP'
332 character*320 controlcard
333 character*3200 controlcard1
334 integer iremd_m_total
336 if(me.eq.king.or..not.out1file)
337 & write (iout,*) "REMD setup"
339 call card_concat(controlcard)
340 call readi(controlcard,"NREP",nrep,3)
341 call readi(controlcard,"NSTEX",nstex,1000)
342 call reada(controlcard,"RETMIN",retmin,10.0d0)
343 call reada(controlcard,"RETMAX",retmax,1000.0d0)
344 mremdsync=(index(controlcard,'SYNC').gt.0)
345 call readi(controlcard,"NSYN",i_sync_step,100)
346 restart1file=(index(controlcard,'REST1FILE').gt.0)
347 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
348 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
349 if(max_cache_traj_use.gt.max_cache_traj)
350 & max_cache_traj_use=max_cache_traj
351 if(me.eq.king.or..not.out1file) then
352 cd if (traj1file) then
353 crc caching is in testing - NTWX is not ignored
354 cd write (iout,*) "NTWX value is ignored"
355 cd write (iout,*) " trajectory is stored to one file by master"
356 cd write (iout,*) " before exchange at NSTEX intervals"
358 write (iout,*) "NREP= ",nrep
359 write (iout,*) "NSTEX= ",nstex
360 write (iout,*) "SYNC= ",mremdsync
361 write (iout,*) "NSYN= ",i_sync_step
362 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
365 if (index(controlcard,'TLIST').gt.0) then
367 call card_concat(controlcard1)
368 read(controlcard1,*) (remd_t(i),i=1,nrep)
369 if(me.eq.king.or..not.out1file)
370 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
373 if (index(controlcard,'MLIST').gt.0) then
375 call card_concat(controlcard1)
376 read(controlcard1,*) (remd_m(i),i=1,nrep)
377 if(me.eq.king.or..not.out1file) then
378 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
381 iremd_m_total=iremd_m_total+remd_m(i)
383 write (iout,*) 'Total number of replicas ',iremd_m_total
388 & "Adaptive (PMF-biased) umbrella sampling will be run"
391 if(me.eq.king.or..not.out1file)
392 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
395 c--------------------------------------------------------------------------
396 subroutine read_MDpar
400 implicit real*8 (a-h,o-z)
402 include 'COMMON.IOUNITS'
403 include 'COMMON.TIME1'
406 include 'COMMON.LANGEVIN'
408 include 'COMMON.LANGEVIN.lang0'
410 include 'COMMON.INTERACT'
411 include 'COMMON.NAMES'
413 include 'COMMON.SETUP'
414 include 'COMMON.CONTROL'
415 include 'COMMON.SPLITELE'
416 include 'COMMON.FFIELD'
418 character*320 controlcard
420 call card_concat(controlcard)
421 call readi(controlcard,"NSTEP",n_timestep,1000000)
422 call readi(controlcard,"NTWE",ntwe,100)
423 call readi(controlcard,"NTWX",ntwx,1000)
424 call reada(controlcard,"DT",d_time,1.0d-1)
425 call reada(controlcard,"DVMAX",dvmax,2.0d1)
426 call reada(controlcard,"DAMAX",damax,1.0d1)
427 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
428 call readi(controlcard,"LANG",lang,0)
429 RESPA = index(controlcard,"RESPA") .gt. 0
430 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
431 ntime_split0=ntime_split
432 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
433 ntime_split0=ntime_split
434 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
435 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
436 rest = index(controlcard,"REST").gt.0
437 tbf = index(controlcard,"TBF").gt.0
438 usampl = index(controlcard,"USAMPL").gt.0
439 scale_umb = index(controlcard,"SCALE_UMB").gt.0
440 adaptive = index(controlcard,"ADAPTIVE").gt.0
441 mdpdb = index(controlcard,"MDPDB").gt.0
442 call reada(controlcard,"T_BATH",t_bath,300.0d0)
443 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
444 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
445 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
446 if (count_reset_moment.eq.0) count_reset_moment=1000000000
447 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
448 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
449 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
450 if (count_reset_vel.eq.0) count_reset_vel=1000000000
451 large = index(controlcard,"LARGE").gt.0
452 print_compon = index(controlcard,"PRINT_COMPON").gt.0
453 rattle = index(controlcard,"RATTLE").gt.0
454 c if performing umbrella sampling, fragments constrained are read from the fragment file
457 write (iout,*) "Umbrella sampling will be run"
458 if (scale_umb.and.adaptive) then
459 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
460 write (iout,*) "Select one of those and re-run the job."
463 if (scale_umb) write (iout,*)
464 &"Umbrella-restraint force constants will be scaled by temperature"
468 if(me.eq.king.or..not.out1file) then
470 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
472 write (iout,'(a)') "The units are:"
473 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
474 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
475 & " acceleration: angstrom/(48.9 fs)**2"
476 write (iout,'(a)') "energy: kcal/mol, temperature: K"
478 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
479 write (iout,'(a60,f10.5,a)')
480 & "Initial time step of numerical integration:",d_time,
482 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
484 write (iout,'(2a,i4,a)')
485 & "A-MTS algorithm used; initial time step for fast-varying",
486 & " short-range forces split into",ntime_split," steps."
487 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
488 & r_cut," lambda",rlamb
490 write (iout,'(2a,f10.5)')
491 & "Maximum acceleration threshold to reduce the time step",
492 & "/increase split number:",damax
493 write (iout,'(2a,f10.5)')
494 & "Maximum predicted energy drift to reduce the timestep",
495 & "/increase split number:",edriftmax
496 write (iout,'(a60,f10.5)')
497 & "Maximum velocity threshold to reduce velocities:",dvmax
498 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
499 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
500 if (rattle) write (iout,'(a60)')
501 & "Rattle algorithm used to constrain the virtual bonds"
505 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
506 call reada(controlcard,"RWAT",rwat,1.4d0)
507 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
508 surfarea=index(controlcard,"SURFAREA").gt.0
509 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
510 if(me.eq.king.or..not.out1file)then
511 write (iout,'(/a,$)') "Langevin dynamics calculation"
514 & " with direct integration of Langevin equations"
515 else if (lang.eq.2) then
516 write (iout,'(a/)') " with TINKER stochasic MD integrator"
517 else if (lang.eq.3) then
518 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
519 else if (lang.eq.4) then
520 write (iout,'(a/)') " in overdamped mode"
522 write (iout,'(//a,i5)')
523 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
526 write (iout,'(a60,f10.5)') "Temperature:",t_bath
527 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
528 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
529 write (iout,'(a60,f10.5)')
530 & "Scaling factor of the friction forces:",scal_fric
531 if (surfarea) write (iout,'(2a,i10,a)')
532 & "Friction coefficients will be scaled by solvent-accessible",
533 & " surface area every",reset_fricmat," steps."
535 c Calculate friction coefficients and bounds of stochastic forces
536 eta=6*pi*cPoise*etawat
537 if(me.eq.king.or..not.out1file)
538 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
540 gamp=scal_fric*(pstok+rwat)*eta
541 stdfp=dsqrt(2*Rb*t_bath/d_time)
543 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
544 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
546 if(me.eq.king.or..not.out1file)then
547 write (iout,'(/2a/)')
548 & "Radii of site types and friction coefficients and std's of",
549 & " stochastic forces of fully exposed sites"
550 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
552 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
553 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
557 if(me.eq.king.or..not.out1file)then
558 write (iout,'(a)') "Berendsen bath calculation"
559 write (iout,'(a60,f10.5)') "Temperature:",t_bath
560 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
562 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
563 & count_reset_moment," steps"
565 & write (iout,'(a,i10,a)')
566 & "Velocities will be reset at random every",count_reset_vel,
570 if(me.eq.king.or..not.out1file)
571 & write (iout,'(a31)') "Microcanonical mode calculation"
573 if(me.eq.king.or..not.out1file)then
574 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
576 write(iout,*) "MD running with constraints."
577 write(iout,*) "Equilibration time ", eq_time, " mtus."
578 write(iout,*) "Constraining ", nfrag," fragments."
579 write(iout,*) "Length of each fragment, weight and q0:"
581 write (iout,*) "Set of restraints #",iset
583 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
584 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
586 write(iout,*) "constraints between ", npair, "fragments."
587 write(iout,*) "constraint pairs, weights and q0:"
589 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
590 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
592 write(iout,*) "angle constraints within ", nfrag_back,
593 & "backbone fragments."
595 write(iout,*) "fragment, weights, q0:"
597 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
598 & ifrag_back(2,i,iset),
599 & wfrag_back(1,i,iset),qin_back(1,i,iset),
600 & wfrag_back(2,i,iset),qin_back(2,i,iset),
601 & wfrag_back(3,i,iset),qin_back(3,i,iset)
604 write(iout,*) "fragment, weights:"
606 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
607 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
608 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
612 iset=mod(kolor,nset)+1
615 if(me.eq.king.or..not.out1file)
616 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
619 c------------------------------------------------------------------------------
622 C Read molecular data.
624 implicit real*8 (a-h,o-z)
630 include 'COMMON.IOUNITS'
633 include 'COMMON.INTERACT'
634 include 'COMMON.LOCAL'
635 include 'COMMON.NAMES'
636 include 'COMMON.CHAIN'
637 include 'COMMON.FFIELD'
638 include 'COMMON.SBRIDGE'
639 include 'COMMON.HEADER'
640 include 'COMMON.CONTROL'
641 include 'COMMON.DBASE'
642 include 'COMMON.THREAD'
643 include 'COMMON.CONTACTS'
644 include 'COMMON.TORCNSTR'
645 include 'COMMON.TIME1'
646 include 'COMMON.BOUNDS'
648 include 'COMMON.SETUP'
649 include 'COMMON.SHIELD'
650 character*4 sequence(maxres)
652 double precision x(maxvar)
653 character*256 pdbfile
654 character*400 weightcard
655 character*80 weightcard_t,ucase
656 dimension itype_pdb(maxres)
657 common /pizda/ itype_pdb
658 logical seq_comp,fail
659 double precision energia(0:n_ene)
660 double precision secprob(3,maxdih_constr)
666 C Read weights of the subsequent energy terms.
667 call card_concat(weightcard)
668 call reada(weightcard,'WLONG',wlong,1.0D0)
669 call reada(weightcard,'WSC',wsc,wlong)
670 call reada(weightcard,'WSCP',wscp,wlong)
671 call reada(weightcard,'WELEC',welec,1.0D0)
672 call reada(weightcard,'WVDWPP',wvdwpp,welec)
673 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
674 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
675 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
676 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
677 call reada(weightcard,'WTURN3',wturn3,1.0D0)
678 call reada(weightcard,'WTURN4',wturn4,1.0D0)
679 call reada(weightcard,'WTURN6',wturn6,1.0D0)
680 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
681 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
682 call reada(weightcard,'WBOND',wbond,1.0D0)
683 call reada(weightcard,'WTOR',wtor,1.0D0)
684 call reada(weightcard,'WTORD',wtor_d,1.0D0)
685 call reada(weightcard,'WANG',wang,1.0D0)
686 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
687 call reada(weightcard,'SCAL14',scal14,0.4D0)
688 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
689 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
690 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
691 call reada(weightcard,'TEMP0',temp0,300.0d0)
692 call reada(weightcard,'WSHIELD',wshield,1.0d0)
693 call reada(weightcard,'WLT',wliptran,0.0D0)
694 call reada(weightcard,'WTUBE',wtube,1.0D0)
695 if (index(weightcard,'SOFT').gt.0) ipot=6
696 C 12/1/95 Added weight for the multi-body term WCORR
697 call reada(weightcard,'WCORRH',wcorr,1.0D0)
698 if (wcorr4.gt.0.0d0) wcorr=wcorr4
718 if(me.eq.king.or..not.out1file)
719 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
720 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
722 10 format (/'Energy-term weights (unscaled):'//
723 & 'WSCC= ',f10.6,' (SC-SC)'/
724 & 'WSCP= ',f10.6,' (SC-p)'/
725 & 'WELEC= ',f10.6,' (p-p electr)'/
726 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
727 & 'WBOND= ',f10.6,' (stretching)'/
728 & 'WANG= ',f10.6,' (bending)'/
729 & 'WSCLOC= ',f10.6,' (SC local)'/
730 & 'WTOR= ',f10.6,' (torsional)'/
731 & 'WTORD= ',f10.6,' (double torsional)'/
732 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
733 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
734 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
735 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
736 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
737 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
738 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
739 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
740 & 'WTURN6= ',f10.6,' (turns, 6th order)')
741 if(me.eq.king.or..not.out1file)then
742 if (wcorr4.gt.0.0d0) then
743 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
744 & 'between contact pairs of peptide groups'
745 write (iout,'(2(a,f5.3/))')
746 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
747 & 'Range of quenching the correlation terms:',2*delt_corr
748 else if (wcorr.gt.0.0d0) then
749 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
750 & 'between contact pairs of peptide groups'
752 write (iout,'(a,f8.3)')
753 & 'Scaling factor of 1,4 SC-p interactions:',scal14
754 write (iout,'(a,f8.3)')
755 & 'General scaling factor of SC-p interactions:',scalscp
757 r0_corr=cutoff_corr-delt_corr
759 aad(i,1)=scalscp*aad(i,1)
760 aad(i,2)=scalscp*aad(i,2)
761 bad(i,1)=scalscp*bad(i,1)
762 bad(i,2)=scalscp*bad(i,2)
765 call rescale_weights(t_bath)
766 if(me.eq.king.or..not.out1file)
767 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
768 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
770 22 format (/'Energy-term weights (scaled):'//
771 & 'WSCC= ',f10.6,' (SC-SC)'/
772 & 'WSCP= ',f10.6,' (SC-p)'/
773 & 'WELEC= ',f10.6,' (p-p electr)'/
774 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
775 & 'WBOND= ',f10.6,' (stretching)'/
776 & 'WANG= ',f10.6,' (bending)'/
777 & 'WSCLOC= ',f10.6,' (SC local)'/
778 & 'WTOR= ',f10.6,' (torsional)'/
779 & 'WTORD= ',f10.6,' (double torsional)'/
780 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
781 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
782 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
783 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
784 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
785 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
786 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
787 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
788 & 'WTURN6= ',f10.6,' (turns, 6th order)')
789 if(me.eq.king.or..not.out1file)
790 & write (iout,*) "Reference temperature for weights calculation:",
792 call reada(weightcard,"D0CM",d0cm,3.78d0)
793 call reada(weightcard,"AKCM",akcm,15.1d0)
794 call reada(weightcard,"AKTH",akth,11.0d0)
795 call reada(weightcard,"AKCT",akct,12.0d0)
796 call reada(weightcard,"V1SS",v1ss,-1.08d0)
797 call reada(weightcard,"V2SS",v2ss,7.61d0)
798 call reada(weightcard,"V3SS",v3ss,13.7d0)
799 call reada(weightcard,"EBR",ebr,-5.50D0)
800 call reada(weightcard,"ATRISS",atriss,0.301D0)
801 call reada(weightcard,"BTRISS",btriss,0.021D0)
802 call reada(weightcard,"CTRISS",ctriss,1.001D0)
803 call reada(weightcard,"DTRISS",dtriss,1.001D0)
804 write (iout,*) "ATRISS=", atriss
805 write (iout,*) "BTRISS=", btriss
806 write (iout,*) "CTRISS=", ctriss
807 write (iout,*) "DTRISS=", dtriss
808 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
810 dyn_ss_mask(i)=.false.
814 dyn_ssbond_ij(i,j)=1.0d300
817 call reada(weightcard,"HT",Ht,0.0D0)
819 ss_depth=ebr/wsc-0.25*eps(1,1)
820 Ht=Ht/wsc-0.25*eps(1,1)
821 akcm=akcm*wstrain/wsc
822 akth=akth*wstrain/wsc
823 akct=akct*wstrain/wsc
824 v1ss=v1ss*wstrain/wsc
825 v2ss=v2ss*wstrain/wsc
826 v3ss=v3ss*wstrain/wsc
828 if (wstrain.ne.0.0) then
829 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
834 write (iout,*) "wshield,", wshield
835 if(me.eq.king.or..not.out1file) then
836 write (iout,*) "Parameters of the SS-bond potential:"
837 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
839 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
840 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
841 write (iout,*)" HT",Ht
842 print *,'indpdb=',indpdb,' pdbref=',pdbref
844 if (indpdb.gt.0 .or. pdbref) then
845 read(inp,'(a)') pdbfile
846 if(me.eq.king.or..not.out1file)
847 & write (iout,'(2a)') 'PDB data will be read from file ',
848 & pdbfile(:ilen(pdbfile))
849 open(ipdbin,file=pdbfile,status='old',err=33)
851 33 write (iout,'(a)') 'Error opening PDB file.'
854 c write (iout,*) 'Begin reading pdb data'
858 c write (iout,*) 'Finished reading pdb data'
860 if(me.eq.king.or..not.out1file)
861 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
862 & ' nstart_sup=',nstart_sup
864 itype_pdb(i)=itype(i)
868 nct=nstart_sup+nsup-1
869 call contact(.false.,ncont_ref,icont_ref,co)
872 if(me.eq.king.or..not.out1file)
873 & write(iout,*)'Adding sidechains'
877 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
880 do while (fail.and.nsi.le.maxsi)
881 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
884 if(fail) write(iout,*)'Adding sidechain failed for res ',
885 & i,' after ',nsi,' trials'
890 if (indpdb.lt.2) then
891 C Read sequence if not taken from the pdb file.
893 c print *,'nres=',nres
894 if (iscode.gt.0) then
895 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
897 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
899 C Convert sequence to numeric code
901 itype(i)=rescode(i,sequence(i),iscode)
903 C Assign initial virtual bond lengths
909 vbld(i+nres)=dsc(iabs(itype(i)))
910 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
911 c write (iout,*) "i",i," itype",itype(i),
912 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
916 c print '(20i4)',(itype(i),i=1,nres)
919 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
921 if (itype(i).eq.ntyp1) then
925 else if (iabs(itype(i+1)).ne.20) then
927 else if (iabs(itype(i)).ne.20) then
934 if(me.eq.king.or..not.out1file)then
935 write (iout,*) "ITEL"
937 write (iout,*) i,itype(i),itel(i)
939 print *,'Call Read_Bridge.'
943 cd print *,'NNT=',NNT,' NCT=',NCT
944 if (itype(1).eq.ntyp1) nnt=2
945 if (itype(nres).eq.ntyp1) nct=nct-1
946 if (indpdb.eq.1 .or. pdbref) then
947 if(me.eq.king.or..not.out1file)
948 & write (iout,'(a,i3)') 'nsup=',nsup
950 if (nsup.le.(nct-nnt+1)) then
951 do i=0,nct-nnt+1-nsup
952 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
958 & 'Error - sequences to be superposed do not match.'
961 do i=0,nsup-(nct-nnt+1)
962 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
964 nstart_sup=nstart_sup+i
970 & 'Error - sequences to be superposed do not match.'
973 if (nsup.eq.0) nsup=nct-nnt
974 if (nstart_sup.eq.0) nstart_sup=nnt
975 if (nstart_seq.eq.0) nstart_seq=nnt
976 if(me.eq.king.or..not.out1file)
977 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
978 & ' nstart_seq=',nstart_seq
979 c AL 7/13/18 Compute quantities to locate chain
980 do i=nstart_sup,nstart_sup+nsup-2
981 call refsys(i,i+1,i+2,prod(1,1,i+2),prod(1,2,i+2),
982 & prod(1,3,i+2),fail)
985 & "Failed to calculate reference system, residue",i
991 C 8/13/98 Set limits to generating the dihedral angles
996 read (inp,*) ndih_constr
997 write (iout,*) "ndih_constr",ndih_constr
998 if (ndih_constr.gt.0) then
1000 C read (inp,*) ftors
1001 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
1004 if(me.eq.king.or..not.out1file)then
1007 & 'There are',ndih_constr,' restraints on gamma angles.'
1009 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
1016 phi0(i)=deg2rad*phi0(i)
1017 drange(i)=deg2rad*drange(i)
1019 C if(me.eq.king.or..not.out1file)
1020 C & write (iout,*) 'FTORS',ftors
1023 phibound(1,ii) = phi0(i)-drange(i)
1024 phibound(2,ii) = phi0(i)+drange(i)
1027 if (me.eq.king .or. .not.out1file) then
1029 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
1031 write (iout,'(a3,i5,2f10.1)')
1032 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1037 else if (ndih_constr.lt.0) then
1039 call card_concat(weightcard)
1040 call reada(weightcard,"PHIHEL",phihel,50.0D0)
1041 call reada(weightcard,"PHIBET",phibet,180.0D0)
1042 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1043 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1044 call reada(weightcard,"WDIHC",wdihc,0.591D0)
1045 write (iout,*) "Weight of dihedral angle restraints",wdihc
1046 read(inp,'(9x,3f7.3)')
1047 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1048 write (iout,*) "The secprob array"
1050 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1054 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
1055 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
1056 ndih_constr=ndih_constr+1
1057 idih_constr(ndih_constr)=i
1060 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1061 sumv=sumv+vpsipred(j,ndih_constr)
1064 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1066 phibound(1,ndih_constr)=phihel*deg2rad
1067 phibound(2,ndih_constr)=phibet*deg2rad
1068 sdihed(1,ndih_constr)=sigmahel*deg2rad
1069 sdihed(2,ndih_constr)=sigmabet*deg2rad
1073 if(me.eq.king.or..not.out1file)then
1076 & 'There are',ndih_constr,
1077 & ' bimodal restraints on gamma angles.'
1079 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1080 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1081 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1082 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1083 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1084 & (vpsipred(j,i),j=1,3)
1091 C first setting the theta boundaries to 0 to pi
1092 C this mean that there is no energy penalty for any angle occuring this can be applied
1093 C for generate random conformation but is not implemented in this way
1096 C thetabound(2,i)=pi
1098 C begin reading theta constrains this is quartic constrains allowing to
1099 C have smooth second derivative
1100 if (with_theta_constr) then
1101 C with_theta_constr is keyword allowing for occurance of theta constrains
1102 read (inp,*) ntheta_constr
1103 C ntheta_constr is the number of theta constrains
1104 if (ntheta_constr.gt.0) then
1105 C read (inp,*) ftors
1106 read (inp,*) (itheta_constr(i),theta_constr0(i),
1107 & theta_drange(i),for_thet_constr(i),
1108 & i=1,ntheta_constr)
1109 C the above code reads from 1 to ntheta_constr
1110 C itheta_constr(i) residue i for which is theta_constr
1111 C theta_constr0 the global minimum value
1112 C theta_drange is range for which there is no energy penalty
1113 C for_thet_constr is the force constant for quartic energy penalty
1115 if(me.eq.king.or..not.out1file)then
1117 & 'There are',ntheta_constr,' constraints on phi angles.'
1118 do i=1,ntheta_constr
1119 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1121 & for_thet_constr(i)
1124 do i=1,ntheta_constr
1125 theta_constr0(i)=deg2rad*theta_constr0(i)
1126 theta_drange(i)=deg2rad*theta_drange(i)
1128 C if(me.eq.king.or..not.out1file)
1129 C & write (iout,*) 'FTORS',ftors
1130 C do i=1,ntheta_constr
1131 C ii = itheta_constr(i)
1132 C thetabound(1,ii) = phi0(i)-drange(i)
1133 C thetabound(2,ii) = phi0(i)+drange(i)
1135 endif ! ntheta_constr.gt.0
1136 endif! with_theta_constr
1138 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1139 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1140 c--- Zscore rms -------
1141 if (nz_start.eq.0) nz_start=nnt
1142 if (nz_end.eq.0 .and. nsup.gt.0) then
1144 else if (nz_end.eq.0) then
1147 if(me.eq.king.or..not.out1file)then
1148 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1149 write (iout,*) 'IZ_SC=',iz_sc
1151 c----------------------
1154 if (.not.pdbref) then
1155 call read_angles(inp,*38)
1157 38 write (iout,'(a)') 'Error reading reference structure.'
1159 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1160 stop 'Error reading reference structure'
1162 39 call chainbuild_extconf
1164 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1171 cref(j,i,kkk)=c(j,i)
1174 call contact(.true.,ncont_ref,icont_ref,co)
1178 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1180 if (constr_dist.gt.0) call read_dist_constr
1181 write (iout,*) "After read_dist_constr nhpb",nhpb
1182 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1184 if(me.eq.king.or..not.out1file)
1185 & write (iout,*) 'Contact order:',co
1187 if(me.eq.king.or..not.out1file)
1188 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1191 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1193 if(me.eq.king.or..not.out1file)
1194 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1195 & icont_ref(1,i),' ',
1196 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1200 if (indpdb.lt.2 .and. modecalc.ne.2 .and. modecalc.ne.4
1201 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1202 & modecalc.ne.10) then
1203 C If input structure hasn't been supplied from the PDB file read or generate
1205 if (iranconf.eq.0 .and. indpdb.eq.0 .and. .not. extconf) then
1206 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1207 & write (iout,'(a)') 'Initial geometry will be read in.'
1209 read(inp,'(8f10.5)',end=36,err=36)
1210 & ((c(l,k),l=1,3),k=1,nres),
1211 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1212 write (iout,*) "Exit READ_CART"
1213 c write (iout,'(8f10.5)')
1214 c & ((c(l,k),l=1,3),k=1,nres),
1215 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1219 dc(j,i)=c(j,i+1)-c(j,i)
1220 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1224 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1226 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1227 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1232 c dc_norm(j,i+nres)=0.0d0
1236 call int_from_cart1(.true.)
1237 write (iout,*) "Finish INT_TO_CART"
1238 c write (iout,*) "DC"
1240 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1241 c & (dc(j,i+nres),j=1,3)
1247 call read_angles(inp,*36)
1248 call chainbuild_extconf
1251 36 write (iout,'(a)') 'Error reading angle file.'
1253 call mpi_finalize( MPI_COMM_WORLD,IERR )
1255 stop 'Error reading angle file.'
1257 else if (extconf) then
1258 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1259 & write (iout,'(a)') 'Extended chain initial geometry.'
1261 theta(i)=90d0*deg2rad
1264 phi(i)=180d0*deg2rad
1267 alph(i)=110d0*deg2rad
1270 omeg(i)=-120d0*deg2rad
1271 if (itype(i).le.0) omeg(i)=-omeg(i)
1273 call chainbuild_extconf
1275 if (indpdb.eq.0) then
1276 if (me.eq.king.or..not.out1file)
1277 & write (iout,'(a)') 'Random-generated initial geometry.'
1279 if (me.eq.king.or..not.out1file)write (iout,'(a)')
1280 &'The geometry of missing residues in PDB file randomly generated.'
1284 if (me.eq.king .or. fg_rank.eq.0 .and. (
1285 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1288 if (indpdb.eq.0) then
1291 call gen_rand_conf(itmp,nrestemp,*30)
1295 do i=nstart_seq,nstart_seq+nsup
1297 c(j,i)=cref(j,i-nstart_seq+nnt,1)
1298 c(j,i+nres)=cref(j,i-nstart_seq+nnt+nres_pdb,1)
1302 & "Coordinates copied from the reference structure"
1304 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
1305 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
1306 & (c(j,nres+ires),j=1,3)
1308 if (nstart_seq.gt.nnt) then
1309 nstartmp=nstart_seq-1
1311 write (iout,*) "Building a random chain from residue",
1312 & nstartmp," down to residue",itmp
1313 call gen_rand_conf(nstartmp,itmp,*30)
1316 itmp=nstart_seq+nsup+1
1317 if (itmp.lt.nres) then
1319 write (iout,*)"Building a random chain from residue",
1320 & itmp," up to residue",nrestemp
1321 call gen_rand_conf(itmp,nrestemp,*30)
1325 30 write (iout,*) 'Failed to generate random conformation',
1326 & ', itrial=',itrial
1327 write (*,*) 'Processor:',me,
1328 & ' Failed to generate random conformation',
1337 write (iout,'(a,i3,a)') 'Processor:',me,
1338 & ' error in generating random conformation.'
1339 write (*,'(a,i3,a)') 'Processor:',me,
1340 & ' error in generating random conformation.'
1343 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1348 & ' error in generating random conformation.'
1352 write (iout,*) "Coordinates after building missing part"
1354 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
1355 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
1356 & (c(j,nres+ires),j=1,3)
1358 call int_from_cart1(.true.)
1360 elseif (modecalc.eq.4) then
1361 read (inp,'(a)') intinname
1362 open (intin,file=intinname,status='old',err=333)
1363 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1364 & write (iout,'(a)') 'intinname',intinname
1365 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1367 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1369 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1371 stop 'Error opening angle file.'
1375 C Generate distance constraints, if the PDB structure is to be regularized.
1376 if (nthread.gt.0) then
1377 call read_threadbase
1380 if (me.eq.king .or. .not. out1file)
1382 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1383 write (iout,'(/a,i3,a)')
1384 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1385 write (iout,'(20i4)') (iss(i),i=1,ns)
1387 write(iout,*)"Running with dynamic disulfide-bond formation"
1389 write (iout,'(/a/)') 'Pre-formed links are:'
1395 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1396 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1402 if (ns.gt.0.and.dyn_ss) then
1406 forcon(i-nss)=forcon(i)
1413 dyn_ss_mask(iss(i))=.true.
1418 c write (iout,*) "DC"
1420 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1421 c & (dc(j,i+nres),j=1,3)
1423 if (i2ndstr.gt.0) call secstrp2dihc
1424 c call geom_to_var(nvar,x)
1425 c call etotal(energia(0))
1426 c call enerprint(energia(0))
1427 c call briefout(0,etot)
1429 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1430 cd write (iout,'(a)') 'Variable list:'
1431 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1433 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1434 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1435 & 'Processor',myrank,': end reading molecular data.'
1440 c--------------------------------------------------------------------------
1441 logical function seq_comp(itypea,itypeb,length)
1443 integer length,itypea(length),itypeb(length)
1446 if (itypea(i).ne.itypeb(i)) then
1454 c-----------------------------------------------------------------------------
1455 subroutine read_bridge
1456 C Read information about disulfide bridges.
1457 implicit real*8 (a-h,o-z)
1458 include 'DIMENSIONS'
1462 include 'COMMON.IOUNITS'
1463 include 'COMMON.GEO'
1464 include 'COMMON.VAR'
1465 include 'COMMON.INTERACT'
1466 include 'COMMON.LOCAL'
1467 include 'COMMON.NAMES'
1468 include 'COMMON.CHAIN'
1469 include 'COMMON.FFIELD'
1470 include 'COMMON.SBRIDGE'
1471 include 'COMMON.HEADER'
1472 include 'COMMON.CONTROL'
1473 include 'COMMON.DBASE'
1474 include 'COMMON.THREAD'
1475 include 'COMMON.TIME1'
1476 include 'COMMON.SETUP'
1477 C Read bridging residues.
1478 read (inp,*) ns,(iss(i),i=1,ns)
1480 if(me.eq.king.or..not.out1file)
1481 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1482 C Check whether the specified bridging residues are cystines.
1484 if (iabs(itype(iss(i))).ne.1) then
1485 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1486 & 'Do you REALLY think that the residue ',
1487 & restyp(itype(iss(i))),i,
1488 & ' can form a disulfide bridge?!!!'
1489 write (*,'(2a,i3,a)')
1490 & 'Do you REALLY think that the residue ',
1491 & restyp(itype(iss(i))),i,
1492 & ' can form a disulfide bridge?!!!'
1494 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1499 C Read preformed bridges.
1501 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1503 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1506 C Check if the residues involved in bridges are in the specified list of
1507 C bridging residues.
1510 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1511 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1512 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1513 & ' contains residues present in other pairs.'
1514 write (*,'(a,i3,a)') 'Disulfide pair',i,
1515 & ' contains residues present in other pairs.'
1517 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1523 if (ihpb(i).eq.iss(j)) goto 10
1525 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1528 if (jhpb(i).eq.iss(j)) goto 20
1530 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1536 ihpb(i)=ihpb(i)+nres
1537 jhpb(i)=jhpb(i)+nres
1543 c----------------------------------------------------------------------------
1544 subroutine read_x(kanal,*)
1545 implicit real*8 (a-h,o-z)
1546 include 'DIMENSIONS'
1547 include 'COMMON.GEO'
1548 include 'COMMON.VAR'
1549 include 'COMMON.CHAIN'
1550 include 'COMMON.IOUNITS'
1551 include 'COMMON.CONTROL'
1552 include 'COMMON.LOCAL'
1553 include 'COMMON.INTERACT'
1554 c Read coordinates from input
1556 read(kanal,'(8f10.5)',end=10,err=10)
1557 & ((c(l,k),l=1,3),k=1,nres),
1558 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1561 c(j,2*nres)=c(j,nres)
1563 call int_from_cart1(.false.)
1566 dc(j,i)=c(j,i+1)-c(j,i)
1567 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1571 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1573 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1574 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1582 c----------------------------------------------------------------------------
1583 subroutine read_threadbase
1584 implicit real*8 (a-h,o-z)
1585 include 'DIMENSIONS'
1586 include 'COMMON.IOUNITS'
1587 include 'COMMON.GEO'
1588 include 'COMMON.VAR'
1589 include 'COMMON.INTERACT'
1590 include 'COMMON.LOCAL'
1591 include 'COMMON.NAMES'
1592 include 'COMMON.CHAIN'
1593 include 'COMMON.FFIELD'
1594 include 'COMMON.SBRIDGE'
1595 include 'COMMON.HEADER'
1596 include 'COMMON.CONTROL'
1597 include 'COMMON.DBASE'
1598 include 'COMMON.THREAD'
1599 include 'COMMON.TIME1'
1600 C Read pattern database for threading.
1601 read (icbase,*) nseq
1603 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1604 & nres_base(2,i),nres_base(3,i)
1605 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1607 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1608 c & nres_base(2,i),nres_base(3,i)
1609 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1613 if (weidis.eq.0.0D0) weidis=0.1D0
1622 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1623 write (iout,'(a,i5)') 'nexcl: ',nexcl
1624 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1627 c------------------------------------------------------------------------------
1628 subroutine setup_var
1629 implicit real*8 (a-h,o-z)
1630 include 'DIMENSIONS'
1631 include 'COMMON.IOUNITS'
1632 include 'COMMON.GEO'
1633 include 'COMMON.VAR'
1634 include 'COMMON.INTERACT'
1635 include 'COMMON.LOCAL'
1636 include 'COMMON.NAMES'
1637 include 'COMMON.CHAIN'
1638 include 'COMMON.FFIELD'
1639 include 'COMMON.SBRIDGE'
1640 include 'COMMON.HEADER'
1641 include 'COMMON.CONTROL'
1642 include 'COMMON.DBASE'
1643 include 'COMMON.THREAD'
1644 include 'COMMON.TIME1'
1645 C Set up variable list.
1650 write (iout,*) "SETUP_VAR ialph"
1652 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1654 ialph(i,1)=nvar+nside
1658 if (indphi.gt.0) then
1660 else if (indback.gt.0) then
1665 write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1668 c----------------------------------------------------------------------------
1669 subroutine gen_dist_constr
1670 C Generate CA distance constraints.
1671 implicit real*8 (a-h,o-z)
1672 include 'DIMENSIONS'
1673 include 'COMMON.IOUNITS'
1674 include 'COMMON.GEO'
1675 include 'COMMON.VAR'
1676 include 'COMMON.INTERACT'
1677 include 'COMMON.LOCAL'
1678 include 'COMMON.NAMES'
1679 include 'COMMON.CHAIN'
1680 include 'COMMON.FFIELD'
1681 include 'COMMON.SBRIDGE'
1682 include 'COMMON.HEADER'
1683 include 'COMMON.CONTROL'
1684 include 'COMMON.DBASE'
1685 include 'COMMON.THREAD'
1686 include 'COMMON.TIME1'
1687 dimension itype_pdb(maxres)
1688 common /pizda/ itype_pdb
1690 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1691 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1692 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1694 do i=nstart_sup,nstart_sup+nsup-1
1695 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1696 cd & ' seq_pdb', restyp(itype_pdb(i))
1697 do j=i+2,nstart_sup+nsup-1
1699 ihpb(nhpb)=i+nstart_seq-nstart_sup
1700 jhpb(nhpb)=j+nstart_seq-nstart_sup
1702 dhpb(nhpb)=dist(i,j)
1705 cd write (iout,'(a)') 'Distance constraints:'
1710 cd if (ii.gt.nres) then
1715 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1716 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1717 cd & dhpb(i),forcon(i)
1721 c----------------------------------------------------------------------------
1723 implicit real*8 (a-h,o-z)
1724 include 'DIMENSIONS'
1725 include 'COMMON.MAP'
1726 include 'COMMON.IOUNITS'
1727 character*3 angid(4) /'THE','PHI','ALP','OME'/
1728 character*80 mapcard,ucase
1730 read (inp,'(a)') mapcard
1731 mapcard=ucase(mapcard)
1732 if (index(mapcard,'PHI').gt.0) then
1734 else if (index(mapcard,'THE').gt.0) then
1736 else if (index(mapcard,'ALP').gt.0) then
1738 else if (index(mapcard,'OME').gt.0) then
1741 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1742 stop 'Error - illegal variable spec in MAP card.'
1744 call readi (mapcard,'RES1',res1(imap),0)
1745 call readi (mapcard,'RES2',res2(imap),0)
1746 if (res1(imap).eq.0) then
1747 res1(imap)=res2(imap)
1748 else if (res2(imap).eq.0) then
1749 res2(imap)=res1(imap)
1751 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1753 & 'Error - illegal definition of variable group in MAP.'
1754 stop 'Error - illegal definition of variable group in MAP.'
1756 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1757 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1758 call readi(mapcard,'NSTEP',nstep(imap),0)
1759 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1761 & 'Illegal boundary and/or step size specification in MAP.'
1762 stop 'Illegal boundary and/or step size specification in MAP.'
1767 c----------------------------------------------------------------------------
1769 implicit real*8 (a-h,o-z)
1770 include 'DIMENSIONS'
1771 include 'COMMON.IOUNITS'
1772 include 'COMMON.GEO'
1773 include 'COMMON.CSA'
1774 include 'COMMON.BANK'
1775 include 'COMMON.CONTROL'
1777 character*620 mcmcard
1778 call card_concat(mcmcard)
1780 call readi(mcmcard,'NCONF',nconf,50)
1781 call readi(mcmcard,'NADD',nadd,0)
1782 call readi(mcmcard,'JSTART',jstart,1)
1783 call readi(mcmcard,'JEND',jend,1)
1784 call readi(mcmcard,'NSTMAX',nstmax,500000)
1785 call readi(mcmcard,'N0',n0,1)
1786 call readi(mcmcard,'N1',n1,6)
1787 call readi(mcmcard,'N2',n2,4)
1788 call readi(mcmcard,'N3',n3,0)
1789 call readi(mcmcard,'N4',n4,0)
1790 call readi(mcmcard,'N5',n5,0)
1791 call readi(mcmcard,'N6',n6,10)
1792 call readi(mcmcard,'N7',n7,0)
1793 call readi(mcmcard,'N8',n8,0)
1794 call readi(mcmcard,'N9',n9,0)
1795 call readi(mcmcard,'N14',n14,0)
1796 call readi(mcmcard,'N15',n15,0)
1797 call readi(mcmcard,'N16',n16,0)
1798 call readi(mcmcard,'N17',n17,0)
1799 call readi(mcmcard,'N18',n18,0)
1801 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1803 call readi(mcmcard,'NDIFF',ndiff,2)
1804 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1805 call readi(mcmcard,'IS1',is1,1)
1806 call readi(mcmcard,'IS2',is2,8)
1807 call readi(mcmcard,'NRAN0',nran0,4)
1808 call readi(mcmcard,'NRAN1',nran1,2)
1809 call readi(mcmcard,'IRR',irr,1)
1810 call readi(mcmcard,'NSEED',nseed,20)
1811 call readi(mcmcard,'NTOTAL',ntotal,10000)
1812 call reada(mcmcard,'CUT1',cut1,2.0d0)
1813 call reada(mcmcard,'CUT2',cut2,5.0d0)
1814 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1815 call readi(mcmcard,'ICMAX',icmax,3)
1816 call readi(mcmcard,'IRESTART',irestart,0)
1817 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1820 call reada(mcmcard,'DELE',dele,20.0d0)
1821 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1822 call readi(mcmcard,'IREF',iref,0)
1823 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1824 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1825 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1826 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1827 write (iout,*) "NCONF_IN",nconf_in
1830 c----------------------------------------------------------------------------
1831 cfmc subroutine mcmfread
1832 cfmc implicit real*8 (a-h,o-z)
1833 cfmc include 'DIMENSIONS'
1834 cfmc include 'COMMON.MCMF'
1835 cfmc include 'COMMON.IOUNITS'
1836 cfmc include 'COMMON.GEO'
1837 cfmc character*80 ucase
1838 cfmc character*620 mcmcard
1839 cfmc call card_concat(mcmcard)
1841 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1842 cfmc write(iout,*)'MAXRANT=',maxrant
1843 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1844 cfmc write(iout,*)'MAXFAM=',maxfam
1845 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1846 cfmc write(iout,*)'NNET1=',nnet1
1847 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1848 cfmc write(iout,*)'NNET2=',nnet2
1849 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1850 cfmc write(iout,*)'NNET3=',nnet3
1851 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1852 cfmc write(iout,*)'ILASTT=',ilastt
1853 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1854 cfmc write(iout,*)'MAXSTR=',maxstr
1855 cfmc maxstr_f=maxstr/maxfam
1856 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1857 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1858 cfmc write(iout,*)'NMCMF=',nmcmf
1859 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1860 cfmc write(iout,*)'IFOCUS=',ifocus
1861 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1862 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1863 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1864 cfmc write(iout,*)'INTPRT=',intprt
1865 cfmc call readi(mcmcard,'IPRT',iprt,100)
1866 cfmc write(iout,*)'IPRT=',iprt
1867 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1868 cfmc write(iout,*)'IMAXTR=',imaxtr
1869 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1870 cfmc write(iout,*)'MAXEVEN=',maxeven
1871 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1872 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1873 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1874 cfmc write(iout,*)'INIMIN=',inimin
1875 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1876 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1877 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1878 cfmc write(iout,*)'NTHREAD=',nthread
1879 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1880 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1881 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1882 cfmc write(iout,*)'MAXPERT=',maxpert
1883 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1884 cfmc write(iout,*)'IRMSD=',irmsd
1885 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1886 cfmc write(iout,*)'DENEMIN=',denemin
1887 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1888 cfmc write(iout,*)'RCUT1S=',rcut1s
1889 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1890 cfmc write(iout,*)'RCUT1E=',rcut1e
1891 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1892 cfmc write(iout,*)'RCUT2S=',rcut2s
1893 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1894 cfmc write(iout,*)'RCUT2E=',rcut2e
1895 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1896 cfmc write(iout,*)'DPERT1=',d_pert1
1897 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1898 cfmc write(iout,*)'DPERT1A=',d_pert1a
1899 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1900 cfmc write(iout,*)'DPERT2=',d_pert2
1901 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1902 cfmc write(iout,*)'DPERT2A=',d_pert2a
1903 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1904 cfmc write(iout,*)'DPERT2B=',d_pert2b
1905 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1906 cfmc write(iout,*)'DPERT2C=',d_pert2c
1907 cfmc d_pert1=deg2rad*d_pert1
1908 cfmc d_pert1a=deg2rad*d_pert1a
1909 cfmc d_pert2=deg2rad*d_pert2
1910 cfmc d_pert2a=deg2rad*d_pert2a
1911 cfmc d_pert2b=deg2rad*d_pert2b
1912 cfmc d_pert2c=deg2rad*d_pert2c
1913 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1914 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1915 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1916 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1917 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1918 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1919 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1920 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1921 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1922 cfmc write(iout,*)'RCUTINI=',rcutini
1923 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1924 cfmc write(iout,*)'GRAT=',grat
1925 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1926 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1930 c----------------------------------------------------------------------------
1932 implicit real*8 (a-h,o-z)
1933 include 'DIMENSIONS'
1934 include 'COMMON.MCM'
1935 include 'COMMON.MCE'
1936 include 'COMMON.IOUNITS'
1938 character*320 mcmcard
1939 call card_concat(mcmcard)
1940 call readi(mcmcard,'MAXACC',maxacc,100)
1941 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1942 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1943 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1944 call readi(mcmcard,'MAXREPM',maxrepm,200)
1945 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1946 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1947 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1948 call reada(mcmcard,'E_UP',e_up,5.0D0)
1949 call reada(mcmcard,'DELTE',delte,0.1D0)
1950 call readi(mcmcard,'NSWEEP',nsweep,5)
1951 call readi(mcmcard,'NSTEPH',nsteph,0)
1952 call readi(mcmcard,'NSTEPC',nstepc,0)
1953 call reada(mcmcard,'TMIN',tmin,298.0D0)
1954 call reada(mcmcard,'TMAX',tmax,298.0D0)
1955 call readi(mcmcard,'NWINDOW',nwindow,0)
1956 call readi(mcmcard,'PRINT_MC',print_mc,0)
1957 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1958 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1959 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1960 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1961 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1962 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1963 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1964 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1965 if (nwindow.gt.0) then
1966 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1968 winlen(i)=winend(i)-winstart(i)+1
1971 if (tmax.lt.tmin) tmax=tmin
1972 if (tmax.eq.tmin) then
1976 if (nstepc.gt.0 .and. nsteph.gt.0) then
1977 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1978 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1980 C Probabilities of different move types
1981 sumpro_type(0)=0.0D0
1982 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1983 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1984 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1985 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1986 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1987 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1988 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1990 print *,'i',i,' sumprotype',sumpro_type(i)
1991 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1992 print *,'i',i,' sumprotype',sumpro_type(i)
1996 c----------------------------------------------------------------------------
1997 subroutine read_minim
1998 implicit real*8 (a-h,o-z)
1999 include 'DIMENSIONS'
2000 include 'COMMON.MINIM'
2001 include 'COMMON.IOUNITS'
2003 character*320 minimcard
2004 call card_concat(minimcard)
2005 call readi(minimcard,'MAXMIN',maxmin,2000)
2006 call readi(minimcard,'MAXFUN',maxfun,5000)
2007 call readi(minimcard,'MINMIN',minmin,maxmin)
2008 call readi(minimcard,'MINFUN',minfun,maxmin)
2009 call reada(minimcard,'TOLF',tolf,1.0D-2)
2010 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
2011 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
2012 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
2013 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
2014 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2015 & 'Options in energy minimization:'
2016 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
2017 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
2018 & 'MinMin:',MinMin,' MinFun:',MinFun,
2019 & ' TolF:',TolF,' RTolF:',RTolF
2022 c----------------------------------------------------------------------------
2023 subroutine read_angles(kanal,*)
2024 implicit real*8 (a-h,o-z)
2025 include 'DIMENSIONS'
2026 include 'COMMON.GEO'
2027 include 'COMMON.VAR'
2028 include 'COMMON.CHAIN'
2029 include 'COMMON.IOUNITS'
2030 include 'COMMON.CONTROL'
2031 c Read angles from input
2033 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2034 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2035 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2036 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2039 c 9/7/01 avoid 180 deg valence angle
2040 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2042 theta(i)=deg2rad*theta(i)
2043 phi(i)=deg2rad*phi(i)
2044 alph(i)=deg2rad*alph(i)
2045 omeg(i)=deg2rad*omeg(i)
2050 c----------------------------------------------------------------------------
2051 subroutine reada(rekord,lancuch,wartosc,default)
2053 character*(*) rekord,lancuch
2054 double precision wartosc,default
2057 iread=index(rekord,lancuch)
2058 if (iread.eq.0) then
2062 iread=iread+ilen(lancuch)+1
2063 read (rekord(iread:),*,err=10,end=10) wartosc
2068 c----------------------------------------------------------------------------
2069 subroutine readi(rekord,lancuch,wartosc,default)
2071 character*(*) rekord,lancuch
2072 integer wartosc,default
2075 iread=index(rekord,lancuch)
2076 if (iread.eq.0) then
2080 iread=iread+ilen(lancuch)+1
2081 read (rekord(iread:),*,err=10,end=10) wartosc
2086 c----------------------------------------------------------------------------
2087 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2090 integer tablica(dim),default
2091 character*(*) rekord,lancuch
2098 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2099 if (iread.eq.0) return
2100 iread=iread+ilen(lancuch)+1
2101 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2104 c----------------------------------------------------------------------------
2105 subroutine multreada(rekord,lancuch,tablica,dim,default)
2108 double precision tablica(dim),default
2109 character*(*) rekord,lancuch
2116 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2117 if (iread.eq.0) return
2118 iread=iread+ilen(lancuch)+1
2119 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2122 c----------------------------------------------------------------------------
2123 subroutine openunits
2124 implicit real*8 (a-h,o-z)
2125 include 'DIMENSIONS'
2128 character*16 form,nodename
2131 include 'COMMON.SETUP'
2132 include 'COMMON.IOUNITS'
2134 include 'COMMON.CONTROL'
2135 integer lenpre,lenpot,ilen,lentmp
2137 character*3 out1file_text,ucase
2142 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2143 call getenv_loc("PREFIX",prefix)
2145 call getenv_loc("POT",pot)
2146 call getenv_loc("DIRTMP",tmpdir)
2147 call getenv_loc("CURDIR",curdir)
2148 call getenv_loc("OUT1FILE",out1file_text)
2149 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2150 out1file_text=ucase(out1file_text)
2151 if (out1file_text(1:1).eq."Y") then
2154 out1file=fg_rank.gt.0
2159 if (lentmp.gt.0) then
2160 write (*,'(80(1h!))')
2161 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2162 write (*,'(80(1h!))')
2163 write (*,*)"All output files will be on node /tmp directory."
2165 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2166 if (me.eq.king) then
2167 write (*,*) "The master node is ",nodename
2168 else if (fg_rank.eq.0) then
2169 write (*,*) "I am the CG slave node ",nodename
2171 write (*,*) "I am the FG slave node ",nodename
2174 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2175 lenpre = lentmp+lenpre+1
2177 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2178 C Get the names and open the input files
2179 #if defined(WINIFL) || defined(WINPGI)
2180 open(1,file=pref_orig(:ilen(pref_orig))//
2181 & '.inp',status='old',readonly,shared)
2182 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2183 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2184 C Get parameter filenames and open the parameter files.
2185 call getenv_loc('BONDPAR',bondname)
2186 open (ibond,file=bondname,status='old',readonly,shared)
2187 call getenv_loc('THETPAR',thetname)
2188 open (ithep,file=thetname,status='old',readonly,shared)
2189 call getenv_loc('ROTPAR',rotname)
2190 open (irotam,file=rotname,status='old',readonly,shared)
2191 call getenv_loc('TORPAR',torname)
2192 open (itorp,file=torname,status='old',readonly,shared)
2193 call getenv_loc('TORDPAR',tordname)
2194 open (itordp,file=tordname,status='old',readonly,shared)
2195 call getenv_loc('FOURIER',fouriername)
2196 open (ifourier,file=fouriername,status='old',readonly,shared)
2197 call getenv_loc('ELEPAR',elename)
2198 open (ielep,file=elename,status='old',readonly,shared)
2199 call getenv_loc('SIDEPAR',sidename)
2200 open (isidep,file=sidename,status='old',readonly,shared)
2201 call getenv_loc('LIPTRANPAR',liptranname)
2202 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2203 #elif (defined CRAY) || (defined AIX)
2204 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2206 c print *,"Processor",myrank," opened file 1"
2207 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2208 c print *,"Processor",myrank," opened file 9"
2209 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2210 C Get parameter filenames and open the parameter files.
2211 call getenv_loc('BONDPAR',bondname)
2212 open (ibond,file=bondname,status='old',action='read')
2213 c print *,"Processor",myrank," opened file IBOND"
2214 call getenv_loc('THETPAR',thetname)
2215 open (ithep,file=thetname,status='old',action='read')
2217 call getenv_loc('THETPARPDB',thetname_pdb)
2218 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2220 c print *,"Processor",myrank," opened file ITHEP"
2221 call getenv_loc('ROTPAR',rotname)
2222 open (irotam,file=rotname,status='old',action='read')
2224 call getenv_loc('ROTPARPDB',rotname_pdb)
2225 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2227 c print *,"Processor",myrank," opened file IROTAM"
2228 call getenv_loc('TORPAR',torname)
2229 open (itorp,file=torname,status='old',action='read')
2230 c print *,"Processor",myrank," opened file ITORP"
2231 call getenv_loc('TORDPAR',tordname)
2232 open (itordp,file=tordname,status='old',action='read')
2233 c print *,"Processor",myrank," opened file ITORDP"
2234 call getenv_loc('SCCORPAR',sccorname)
2235 open (isccor,file=sccorname,status='old',action='read')
2236 c print *,"Processor",myrank," opened file ISCCOR"
2237 call getenv_loc('FOURIER',fouriername)
2238 open (ifourier,file=fouriername,status='old',action='read')
2239 c print *,"Processor",myrank," opened file IFOURIER"
2240 call getenv_loc('ELEPAR',elename)
2241 open (ielep,file=elename,status='old',action='read')
2242 c print *,"Processor",myrank," opened file IELEP"
2243 call getenv_loc('SIDEPAR',sidename)
2244 open (isidep,file=sidename,status='old',action='read')
2245 call getenv_loc('LIPTRANPAR',liptranname)
2246 open (iliptranpar,file=liptranname,status='old',action='read')
2247 c print *,"Processor",myrank," opened file ISIDEP"
2248 c print *,"Processor",myrank," opened parameter files"
2250 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2251 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2252 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2253 C Get parameter filenames and open the parameter files.
2254 call getenv_loc('BONDPAR',bondname)
2255 open (ibond,file=bondname,status='old')
2256 call getenv_loc('THETPAR',thetname)
2257 open (ithep,file=thetname,status='old')
2259 call getenv_loc('THETPARPDB',thetname_pdb)
2260 open (ithep_pdb,file=thetname_pdb,status='old')
2262 call getenv_loc('ROTPAR',rotname)
2263 open (irotam,file=rotname,status='old')
2265 call getenv_loc('ROTPARPDB',rotname_pdb)
2266 open (irotam_pdb,file=rotname_pdb,status='old')
2268 call getenv_loc('TORPAR',torname)
2269 open (itorp,file=torname,status='old')
2270 call getenv_loc('TORDPAR',tordname)
2271 open (itordp,file=tordname,status='old')
2272 call getenv_loc('SCCORPAR',sccorname)
2273 open (isccor,file=sccorname,status='old')
2274 call getenv_loc('FOURIER',fouriername)
2275 open (ifourier,file=fouriername,status='old')
2276 call getenv_loc('ELEPAR',elename)
2277 open (ielep,file=elename,status='old')
2278 call getenv_loc('SIDEPAR',sidename)
2279 open (isidep,file=sidename,status='old')
2280 call getenv_loc('LIPTRANPAR',liptranname)
2281 open (iliptranpar,file=liptranname,status='old')
2283 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2285 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2286 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2287 C Get parameter filenames and open the parameter files.
2288 call getenv_loc('BONDPAR',bondname)
2289 open (ibond,file=bondname,status='old',readonly)
2290 call getenv_loc('THETPAR',thetname)
2291 open (ithep,file=thetname,status='old',readonly)
2292 call getenv_loc('ROTPAR',rotname)
2293 open (irotam,file=rotname,status='old',readonly)
2294 call getenv_loc('TORPAR',torname)
2295 open (itorp,file=torname,status='old',readonly)
2296 call getenv_loc('TORDPAR',tordname)
2297 open (itordp,file=tordname,status='old',readonly)
2298 call getenv_loc('SCCORPAR',sccorname)
2299 open (isccor,file=sccorname,status='old',readonly)
2301 call getenv_loc('THETPARPDB',thetname_pdb)
2302 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2304 call getenv_loc('FOURIER',fouriername)
2305 open (ifourier,file=fouriername,status='old',readonly)
2306 call getenv_loc('ELEPAR',elename)
2307 open (ielep,file=elename,status='old',readonly)
2308 call getenv_loc('SIDEPAR',sidename)
2309 open (isidep,file=sidename,status='old',readonly)
2310 call getenv_loc('LIPTRANPAR',liptranname)
2311 open (iliptranpar,file=liptranname,status='old',action='read')
2313 call getenv_loc('ROTPARPDB',rotname_pdb)
2314 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2318 call getenv_loc('TUBEPAR',tubename)
2319 #if defined(WINIFL) || defined(WINPGI)
2320 open (itube,file=tubename,status='old',readonly,shared)
2321 #elif (defined CRAY) || (defined AIX)
2322 open (itube,file=tubename,status='old',action='read')
2324 open (itube,file=tubename,status='old')
2326 open (itube,file=tubename,status='old',readonly)
2331 C 8/9/01 In the newest version SCp interaction constants are read from a file
2332 C Use -DOLDSCP to use hard-coded constants instead.
2334 call getenv_loc('SCPPAR',scpname)
2335 #if defined(WINIFL) || defined(WINPGI)
2336 open (iscpp,file=scpname,status='old',readonly,shared)
2337 #elif (defined CRAY) || (defined AIX)
2338 open (iscpp,file=scpname,status='old',action='read')
2340 open (iscpp,file=scpname,status='old')
2342 open (iscpp,file=scpname,status='old',readonly)
2345 call getenv_loc('PATTERN',patname)
2346 #if defined(WINIFL) || defined(WINPGI)
2347 open (icbase,file=patname,status='old',readonly,shared)
2348 #elif (defined CRAY) || (defined AIX)
2349 open (icbase,file=patname,status='old',action='read')
2351 open (icbase,file=patname,status='old')
2353 open (icbase,file=patname,status='old',readonly)
2356 C Open output file only for CG processes
2357 c print *,"Processor",myrank," fg_rank",fg_rank
2358 if (fg_rank.eq.0) then
2360 if (nodes.eq.1) then
2363 npos = dlog10(dfloat(nodes-1))+1
2365 if (npos.lt.3) npos=3
2366 write (liczba,'(i1)') npos
2367 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2369 write (liczba,form) me
2370 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2371 & liczba(:ilen(liczba))
2372 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2374 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2376 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2377 & liczba(:ilen(liczba))//'.mol2'
2378 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2379 & liczba(:ilen(liczba))//'.stat'
2381 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2382 & //liczba(:ilen(liczba))//'.stat')
2383 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2386 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2387 & liczba(:ilen(liczba))//'.const'
2392 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2393 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2394 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2395 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2396 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2398 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2400 rest2name=prefix(:ilen(prefix))//'.rst'
2402 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2405 #if defined(AIX) || defined(PGI) || defined(CRAY)
2406 if (me.eq.king .or. .not. out1file) then
2407 open(iout,file=outname,status='unknown')
2409 open(iout,file="/dev/null",status="unknown")
2412 if (fg_rank.gt.0) then
2413 write (liczba,'(i3.3)') myrank/nfgtasks
2414 write (ll,'(bz,i3.3)') fg_rank
2415 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2420 open(igeom,file=intname,status='unknown',position='append')
2421 open(ipdb,file=pdbname,status='unknown')
2422 open(imol2,file=mol2name,status='unknown')
2423 open(istat,file=statname,status='unknown',position='append')
2425 c1out open(iout,file=outname,status='unknown')
2428 if (me.eq.king .or. .not.out1file)
2429 & open(iout,file=outname,status='unknown')
2431 if (fg_rank.gt.0) then
2432 write (liczba,'(i3.3)') myrank/nfgtasks
2433 write (ll,'(bz,i3.3)') fg_rank
2434 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2439 open(igeom,file=intname,status='unknown',access='append')
2440 open(ipdb,file=pdbname,status='unknown')
2441 open(imol2,file=mol2name,status='unknown')
2442 open(istat,file=statname,status='unknown',access='append')
2444 c1out open(iout,file=outname,status='unknown')
2447 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2448 csa_seed=prefix(:lenpre)//'.CSA.seed'
2449 csa_history=prefix(:lenpre)//'.CSA.history'
2450 csa_bank=prefix(:lenpre)//'.CSA.bank'
2451 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2452 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2453 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2454 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2455 csa_int=prefix(:lenpre)//'.int'
2456 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2457 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2458 csa_in=prefix(:lenpre)//'.CSA.in'
2459 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2462 write (iout,'(80(1h-))')
2463 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2464 write (iout,'(80(1h-))')
2465 write (iout,*) "Input file : ",
2466 & pref_orig(:ilen(pref_orig))//'.inp'
2467 write (iout,*) "Output file : ",
2468 & outname(:ilen(outname))
2470 write (iout,*) "Sidechain potential file : ",
2471 & sidename(:ilen(sidename))
2473 write (iout,*) "SCp potential file : ",
2474 & scpname(:ilen(scpname))
2476 write (iout,*) "Electrostatic potential file : ",
2477 & elename(:ilen(elename))
2478 write (iout,*) "Cumulant coefficient file : ",
2479 & fouriername(:ilen(fouriername))
2480 write (iout,*) "Torsional parameter file : ",
2481 & torname(:ilen(torname))
2482 write (iout,*) "Double torsional parameter file : ",
2483 & tordname(:ilen(tordname))
2484 write (iout,*) "SCCOR parameter file : ",
2485 & sccorname(:ilen(sccorname))
2486 write (iout,*) "Bond & inertia constant file : ",
2487 & bondname(:ilen(bondname))
2488 write (iout,*) "Bending-torsion parameter file : ",
2489 & thetname(:ilen(thetname))
2490 write (iout,*) "Rotamer parameter file : ",
2491 & rotname(:ilen(rotname))
2492 write (iout,*) "Thetpdb parameter file : ",
2493 & thetname_pdb(:ilen(thetname_pdb))
2494 write (iout,*) "Threading database : ",
2495 & patname(:ilen(patname))
2497 &write (iout,*)" DIRTMP : ",
2499 write (iout,'(80(1h-))')
2503 c----------------------------------------------------------------------------
2504 subroutine card_concat(card)
2505 implicit real*8 (a-h,o-z)
2506 include 'DIMENSIONS'
2507 include 'COMMON.IOUNITS'
2509 character*80 karta,ucase
2511 read (inp,'(a)') karta
2514 do while (karta(80:80).eq.'&')
2515 card=card(:ilen(card)+1)//karta(:79)
2516 read (inp,'(a)') karta
2519 card=card(:ilen(card)+1)//karta
2522 c----------------------------------------------------------------------------------
2524 implicit real*8 (a-h,o-z)
2525 include 'DIMENSIONS'
2526 include 'COMMON.CHAIN'
2527 include 'COMMON.IOUNITS'
2529 open(irest2,file=rest2name,status='unknown')
2530 read(irest2,*) totT,EK,potE,totE,t_bath
2533 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2536 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2539 read (irest2,*) iset
2544 c---------------------------------------------------------------------------------
2545 subroutine read_fragments
2546 implicit real*8 (a-h,o-z)
2547 include 'DIMENSIONS'
2551 include 'COMMON.SETUP'
2552 include 'COMMON.CHAIN'
2553 include 'COMMON.IOUNITS'
2555 include 'COMMON.CONTROL'
2556 read(inp,*) nset,nfrag,npair,nfrag_back
2557 loc_qlike=(nfrag_back.lt.0)
2558 nfrag_back=iabs(nfrag_back)
2559 if(me.eq.king.or..not.out1file)
2560 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2561 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2563 read(inp,*) mset(iset)
2565 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2567 if(me.eq.king.or..not.out1file)
2568 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2569 & ifrag(2,i,iset), qinfrag(i,iset)
2572 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2574 if(me.eq.king.or..not.out1file)
2575 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2576 & ipair(2,i,iset), qinpair(i,iset)
2580 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2581 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2582 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2583 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2584 if(me.eq.king.or..not.out1file)
2585 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2586 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2587 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2588 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2592 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2593 & wfrag_back(3,i,iset),
2594 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2595 if(me.eq.king.or..not.out1file)
2596 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2597 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2603 C---------------------------------------------------------------------------
2604 subroutine read_afminp
2605 implicit real*8 (a-h,o-z)
2606 include 'DIMENSIONS'
2610 include 'COMMON.SETUP'
2611 include 'COMMON.CONTROL'
2612 include 'COMMON.CHAIN'
2613 include 'COMMON.IOUNITS'
2614 include 'COMMON.SBRIDGE'
2615 character*320 afmcard
2617 call card_concat(afmcard)
2618 call readi(afmcard,"BEG",afmbeg,0)
2619 call readi(afmcard,"END",afmend,0)
2620 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2621 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2622 print *,'FORCE=' ,forceAFMconst
2623 CCCC NOW PROPERTIES FOR AFM
2626 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2628 distafminit=dsqrt(distafminit)
2629 print *,'initdist',distafminit
2633 c-------------------------------------------------------------------------------
2634 subroutine read_dist_constr
2635 implicit real*8 (a-h,o-z)
2636 include 'DIMENSIONS'
2640 include 'COMMON.SETUP'
2641 include 'COMMON.CONTROL'
2642 include 'COMMON.CHAIN'
2643 include 'COMMON.IOUNITS'
2644 include 'COMMON.SBRIDGE'
2645 integer ifrag_(2,100),ipair_(2,100)
2646 double precision wfrag_(100),wpair_(100)
2647 character*500 controlcard
2649 c print *, "WCHODZE"
2650 write (iout,*) "Calling read_dist_constr"
2651 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2653 call card_concat(controlcard)
2654 call readi(controlcard,"NFRAG",nfrag_,0)
2655 call readi(controlcard,"NPAIR",npair_,0)
2656 call readi(controlcard,"NDIST",ndist_,0)
2657 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2658 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2659 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2660 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2661 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2662 normalize = index(controlcard,"NORMALIZE").gt.0
2663 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2664 c write (iout,*) "IFRAG"
2666 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2668 c write (iout,*) "IPAIR"
2670 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2674 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2675 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2676 & ifrag_(2,i)=nstart_sup+nsup-1
2677 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2679 if (wfrag_(i).gt.0.0d0) then
2680 do j=ifrag_(1,i),ifrag_(2,i)-1
2681 do k=j+1,ifrag_(2,i)
2682 c write (iout,*) "j",j," k",k
2684 if (constr_dist.eq.1) then
2689 forcon(nhpb)=wfrag_(i)
2690 else if (constr_dist.eq.2) then
2691 if (ddjk.le.dist_cut) then
2696 forcon(nhpb)=wfrag_(i)
2703 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2706 if (.not.out1file .or. me.eq.king)
2707 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2708 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2710 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2711 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2718 if (wpair_(i).gt.0.0d0) then
2726 do j=ifrag_(1,ii),ifrag_(2,ii)
2727 do k=ifrag_(1,jj),ifrag_(2,jj)
2731 forcon(nhpb)=wpair_(i)
2732 dhpb(nhpb)=dist(j,k)
2734 if (.not.out1file .or. me.eq.king)
2735 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2736 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2738 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2739 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2746 write (iout,*) "Distance restraints as read from input"
2748 if (constr_dist.eq.11) then
2749 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2750 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2751 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2752 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2755 if (.not.out1file .or. me.eq.king)
2756 & write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
2757 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2758 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
2760 write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
2761 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2762 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
2764 if (ibecarb(i).gt.0) then
2765 ihpb(i)=ihpb(i)+nres
2766 jhpb(i)=jhpb(i)+nres
2770 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2771 & ibecarb(i),forcon(nhpb+1)
2772 if (forcon(nhpb+1).gt.0.0d0) then
2774 if (ibecarb(i).gt.0) then
2775 ihpb(i)=ihpb(i)+nres
2776 jhpb(i)=jhpb(i)+nres
2778 if (dhpb(nhpb).eq.0.0d0)
2779 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2782 if (.not.out1file .or. me.eq.king)
2783 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2784 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2786 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2787 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2790 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2791 C if (forcon(nhpb+1).gt.0.0d0) then
2793 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2795 if (constr_dist.eq.11 .and. normalize) then
2796 fordepthmax=fordepth(1)
2798 if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i)
2801 fordepth(i)=fordepth(i)/fordepthmax
2803 write (iout,'(/a/4a5,2a8,2a10)')
2804 & "Normalized Lorenzian-like distance restraints",
2805 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V"
2807 write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i),
2808 & dhpb(i),dhpb1(i),forcon(i),fordepth(i)
2815 c-------------------------------------------------------------------------------
2817 subroutine flush(iu)
2822 subroutine flush(iu)
2827 c------------------------------------------------------------------------------
2828 subroutine copy_to_tmp(source)
2829 include "DIMENSIONS"
2830 include "COMMON.IOUNITS"
2831 character*(*) source
2832 character* 256 tmpfile
2836 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2837 inquire(file=tmpfile,exist=ex)
2839 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2840 & " to temporary directory..."
2841 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2842 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2846 c------------------------------------------------------------------------------
2847 subroutine move_from_tmp(source)
2848 include "DIMENSIONS"
2849 include "COMMON.IOUNITS"
2850 character*(*) source
2853 write (*,*) "Moving ",source(:ilen(source)),
2854 & " from temporary directory to working directory"
2855 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2856 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2859 c------------------------------------------------------------------------------
2860 subroutine random_init(seed)
2862 C Initialize random number generator
2864 implicit real*8 (a-h,o-z)
2865 include 'DIMENSIONS'
2868 logical OKRandom, prng_restart
2870 integer iseed_array(4)
2872 include 'COMMON.IOUNITS'
2873 include 'COMMON.TIME1'
2874 include 'COMMON.THREAD'
2875 include 'COMMON.SBRIDGE'
2876 include 'COMMON.CONTROL'
2877 include 'COMMON.MCM'
2878 include 'COMMON.MAP'
2879 include 'COMMON.HEADER'
2880 include 'COMMON.CSA'
2881 include 'COMMON.CHAIN'
2882 include 'COMMON.MUCA'
2884 include 'COMMON.FFIELD'
2885 include 'COMMON.SETUP'
2886 iseed=-dint(dabs(seed))
2887 if (iseed.eq.0) then
2888 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2889 & 'Random seed undefined. The program will stop.'
2890 write (*,'(/80(1h*)/20x,a/80(1h*))')
2891 & 'Random seed undefined. The program will stop.'
2893 call mpi_finalize(mpi_comm_world,ierr)
2895 stop 'Bad random seed.'
2898 if (fg_rank.eq.0) then
2902 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2903 OKRandom = prng_restart(me,iseed)
2906 tmp=65536.0d0**(4-i)
2907 iseed_array(i) = dint(seed/tmp)
2908 seed=seed-iseed_array(i)*tmp
2911 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2912 & (iseed_array(i),i=1,4)
2913 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2914 & (iseed_array(i),i=1,4)
2915 OKRandom = prng_restart(me,iseed_array)
2918 c r1 = prng_next(me)
2919 r1=ran_number(0.0D0,1.0D0)
2921 & write (iout,*) 'ran_num',r1
2922 if (r1.lt.0.0d0) OKRandom=.false.
2924 if (.not.OKRandom) then
2925 write (iout,*) 'PRNG IS NOT WORKING!!!'
2926 print *,'PRNG IS NOT WORKING!!!'
2929 call mpi_abort(mpi_comm_world,error_msg,ierr)
2932 write (iout,*) 'too many processors for parallel prng'
2933 write (*,*) 'too many processors for parallel prng'
2941 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)